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ABSTRACT 


A numerical technique is formulated, in a computer program U2DIIF, for the 
solution of flow over an airfoil executing an arbitrary unsteady motion in an inviscid 
and incompressible medium. The technique extends the well known Panel Methods for 
steady flow into solving a non-linear unsteady flow problem arising from the 
continuous vortex shedding into the trailing wake due to the unsteady motion of the 
airfoil. Numerous case-runs are presented to verify U2DIIF computer code against 
other theoretical and/or numerical methods as well as in cases where limited 
experimental data are obtainable in literatures. These case-runs include airfoils 
undergoing a step change or a modified ramp change of angle-of-attack, airfoils 
executing harmonic oscillation in pitching and plunging motions and _ airfoils 


penetrating a sharp edge gust. 


THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may 
not have been exercised for all cases of interest. While every effort has been made, 
within the time available, to ensure that the programs are free of computational and 
logic errors, they cannot be considered validated. Any application of these programs 


without additional verification is at the risk of the user. 
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I. INTRODUCTION 


A. GENERAL 

In this thesis, a numerical method 1s formulated and coded in a FORTRAN 
computer program, codename U2DIIF (Unsteady 2-Dimensional  Inviscid 
Incompressible Flow), to solve for the flow over an airfoil which is executing an 


unsteady time-dependent motion in an inviscid, incompressible medium. 


B. APPROACH 

The basic approach to this problem is the extension of a very general and 
powerful technique, called Pane! Methods, developed by Hess & Smith [Ref. 1] for 
steady potential flow problems, to include the unsteady motion of the airfoil that 1s 
continuously shedding vorticity into the trailing wake. This vortex shedding process 
creates the non-linearity effects of the problem in that the wake vortices influence the 
flow over the airfoil which in turn alters the vortex shedding as the airfoil proceeds in 
time. [t is this very non-linearity of unsteady flow that distinguishes itself from the well 
known steady Panel Methods solution where the mathematical formulation of the 
problem results in a set of N linear equations in N unknowns which are solved easily 
with the standard Gaussian elimination algorithm. 

The unsteady flow problem is, however, deprived of this relatively easv solution 
technique. Instead, an iterative type of solution is needed for this non-linear problem. 
The correct mathematical model must therefore be formulated to describe the vortex 
shedding process that provides the mechanism for the iteration to proceed towards a 
converged set of solution in each time step. 

It is the objective of this thesis to develop a numerical computer program that 
performs this non-linear potential flow calculation which proceeds step by step in time. 
At each time step, a complete set of potential flow solutions, inclusive of the airfoil 
Oressure Uistribution. force and moment coetlicients. and the trailing vortex wake 


pattern (strengths and positions of shed vortices), is obtained. 


C. SCOPE 
The Panel Method of Hess & Smith, which utilises both the distributed sources 


and vorticities as panel singularities, for steady flow solution is described in Chapter II. 
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Chapter III formulates the mathematical model for the unsteady flow problem 
and its solution procedures, highlighting the essential features in solving the non-linear 
problem of unsteady flow. 

Chapter IV describes the computer program U2DIIF, its essential capabilities, 
limitations and the necessary input set-up for tvpical case-runs. 

The results of some of the case-runs are presented in Chapter V. They are 
compared with other theoretical and/or numerical methods as well as in cases where 
limited experimental data are obtainable in the literature. These case-runs include 
airfoils undergoing a step change or a modified ramp change in angle-of-attack, airfoils 
executing harmonic oscillation in both pitching and plunging motions and airfoils 
penetrating a sharp edge gust. . 

In the concluding remarks of Chapter VI, the future development and application 
potential of this numerical method to other studies of unsteady 2-dimensional inviscid 


incompressible flow are mentioned. 


I 


Il. STEADY FLOW PROBLEM FORMULATION 


A. PRAME OF REFERENCE 

Consider a 2-dimensional airfoul in motion with constant linear velocity — Vg as 
shown in Figure 2.1. Using an (x,y) coordinate system fixed on the airfoil, where the x- 
axis coincides with the chord line originating from the leading edge towards the trailing 
edge of the airfoil, the flow in this frame of reference is steady. That is to sav, the fluid 
velocity and pressure in the flow field depend only on the spatial coordinates (x,y) and 
not on time. The airfoil then appears to be submerged in an onset flow whose velocity 


is Vs, and making an angle of attack, a@, with the x-axis (see Figure 2.1). 


B. STEADY FLOW PANEL METHODS 
1. Definition of Nodes and Panels 

The airfoil surface 1s divided into (n) straight-line segments. called panels, by 
(n+1) arbitrary chosen points, called nodes, distributed over the airfoil contour as 
shown in Figure 2.2. The panel numbering sequence starts with panel 1 on the lower 
surface at the airfoil trailing edge and proceeds clockwise around the airfoil contour so 
that the last panel (panel n) ends on the upper surface, also at the airfoil trailing edge. 

Notice that this numbering sequence dictates that the airfoil body always lies 
on the right hand side of the i" panel as one proceeds from the i" node to the (i+ 1)" 
node. Also the 1%* and the (n+1) nodes coincide at the trailing edge. It therefore 
facilitates, as shown in Figure 2.2, the common definition of unit normal vector n, and 
the unit tangential vector t, for all panels, i.e., n. 1s directed outward from the body into 
the flow and t. is directed from the i™ node to the (i+ 1)" node. 

2. Distribution of Singularities 

Figure 2.2 also indicates that a uniform source distribution q; and a umform 
vorticity distribution 7 are placed on the i" panel. The source strength q; varies from 
panel to panei whereas the vorticity strength ¥ remains the same for all panels. This 
particular choice of singularity distributions is one of the many tvpes of singularity 
combinations (it happened to be the pioneering one though) ever used in a wide variety 
of the so called Panel Methods. The success of representing the flow past an arbitrary 
shaped airfoil by surface singularity distributions lies in the fact that these singularity 


distributions automatically satisfy Laplace’s equation, the governing flow equation for 





Figure 2.1 Frame of Reference for Steady Flow. 
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Figure 2.2. Panel Methods Representation for Steadv Flow. 
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inviscid incompressible flow, and the boundary condition at the far field (9). In 
addition, the superposition principle applies to any linear homogeneous second order 
partial differential equation such as Laplace’s equation. Therefore one can build up an 
overall complicated flow field by the combination of simple flows if the appropriate 
boundary conditions on the airfoil can be satisfied accuratelv. For our case the overail 
flow field (represented by the velocity potential ®) can be built up by three simple 


flows, 

P= 9 +9, + 9, feqmez a 
where @,Q is the potential of the onset flow, 

Doo = Veo (X cosa + y sing) (ccm 2) 


@. is the velocity potential of a source distribution of strength q(s) per unit length, 





42] q(s) 


In r ds (eqn 7.3) 


@, is the velocity potential of a vorticity distribution of strength y(s) per unit length. 


o= — | m 6 ds (eqn 2.4) 





The integrals in Equations 2.3 and 2.4 are performed along the surface 
contour s and (1,8) are polar coordinates of any field point (x,y) measured from the 
airfoil surface at an arbitrary point as shown in Figure 2.3. The difficult task of 
evaluating these integrals has been greatly simplified by our singularity distributions 
postulated to represent the tlow over the airfoil: that is, instead of integrating over the 
entire airfoil contour, we integrate on each panei along a straight line where q, and * 


are constant, then sum up the effects of all panels. Equation 2.1 therefore becomes, 


n 
@M = Voo (x cosa + ysina) + JS nanet j (sh ier = = Q@ ] ds (eqn 2.5) 


ee | 
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(s) 
GY Vo (x cosa + y sind) =o f—Inrds — f oe 


*» 





6 ds 


ad 


* 


Fizuré.2.3 Potéenual Evaluationst a Field Point. 


It can be seen from Equation 2.5 that ® is completely defined if the (n+ 1) 
unknowns, oF (j= 1,2,....,.n) and 7, can be calculated using a numerical technique yet to 
be described. Once the potential @ is solved, the velocity can be evaluated by taking 
grad @. At this point we introduce a definition of disturbance potential, @, as the sum 


of potential due to both the source and vorticity distribution, 
Oo + OL | (eqn 2.6) 
‘Equation 2 Wetheresere reads, 
D = Po + g : (eqn 2.7) 


The total velocity vector is thus, 


eri = V® 
Van + VO 
=Vo + Vo (eqn 2.8) 


The pressure can be obtained from Bernoulli's Equation, 


me V 
Cc, = ——S = 1 - (el? (eqn 2.9) 
P APV en Voo 





Notice that Figure 2.3 indicates that the field point lies off the airfoil surface, 
however, we are interested in field points that are on the airfoil surface. In the case of 


steady flow, the expressions for V and C » are the same for field points lying on or 


total 
off the airfoil surface. It is nevertheless not the same in unsteady flow, as will be seen 


imeemapter li fpinethatv must include the rigid body motion of the airfoil when one 


total 
evaiuates field points on the airfoil surface. 
3. Boundary Conditions 
The boundary conditions to be satisfied include the flow tangency conditions 
and the Kutta Condition. The flow tangency conditions are satisfied at the exterior mid 
points, called control points, of all panels by taking the resultant velocity at each 


control point to have only (V'), but, 


Z| 


(V9). = 0, i=1,2....,n (eqn 2.10) 


where a and C are the tangential and normal components of the total velocity at 
the control point of the i** panel due to the free stream and the velocities induced by 
the source and vorticity distributions on all the panels. j (j= 1,2.,....,n). 

The Kutta condition postulates that the pressures on the upper and lower 
panels at the trailing edge be equal in order that the flow leaves the trailing edge 
smoothly. By using Bernoulli's equation for steady potential flow, this pressure 
equilibrium condition implies that the tangential velocities in the downstream direction 
at the 18¢ and the n™ panel control points must be equal. This fact is certainly 
consistent with the knowledge that when steady flow is established, the total circulation 
over the airfoil does not change if the tangential velocities are the same at the trailing 


edge panels. 
(Vv), a = (Vou (eqn 2.158) 


If one could explicitly express Equations 2.10 and 2.11 in terms of the unknowns 
q, (j= 1,2,....n) and y, the task is then reduced to solving a linear system of (n+1) 


simultaneous equations for the (n+ 1) unknowns. 


C. INFLUENCE COBPFICIENGS 
1. The Concept of Influence Coefficients 

The numerical technique employed in Panel Methods to manipulate 
equations 2.10 and 2.11 into an algebraic system of linear simultaneous equations 
involves the important concept of influence coefficients. An influence coefficient 1s 
defined as the velocity induced at a field point by a unit strength singularity (be it a 
point singularity or a distributed singularity) placed anywhere within the flow field. In 
this case, it is the unit strength singularity distribution on one panel. Recall that 
equations 2.10 and 2.11 simply require the computation of the normal and tangential 
velocity components at all the panei control points. The normal coniponents of 
velocities are essential in satisfying flow tangency conditions while the tangential 
components of velocities are necessary for satisfying the Kutta condition as well as 
computing the pressure distribution. The procedure is thus to compute, at the om panel 
control point, the velocity components induced by the source and _ vorticity 


distributions on all the panels, j (j=1,2,....,.n), including the i“ panel itself. Summation 
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of all the induced velocities. separately for the normal and tangential components, 
together with the free stream velocity components produces ail the required (v"). and 
cn), i= 1,2,....,n. 
2. Notation for Influence Coefficient 
We shall adopt a consistent set of notation for the influence coefficients used 
throughout this documentation. It is so designated to permit easy recognition in that 
each influence coefficient contains all the associated information one needs. An 


influence coefficient is denoted with a superscript and two subsript as follows: 
78 
A Pq 


where ¥ denotes the type of singularity involved, we shall arbitrarily use A, B and C for 
the uniformly distributed source, uniformiv distributed vorticity and point vortex 
respectively. The superscript s is an indicator telling which component the induced 
velocity is. The first subscript p identifies the Meld point where the induced velocitv is 
evaluated. The second subscript q denotes the particular singularity contributing to the 
induced velocity. 

We thus define, for the steady flow problem, the foilowing influence 
coefficients : 


e A®*..: normal velocity component induced at the i panel control point by unit 
strength source distribution on the j panel. 


° .: : tangential velocity component mood at the i panel control point by 
unit Strength source distribution on the j@ ip eho 


e B ..: normal velocity component induced at the i panel control point by unit 
strength vorticity distribution on the j“ panel. 


° Br. ; tangential velocity component nae at the i panel control point by 
unit strength vorticity distribution on the j" eee omel: 


3. Computation of Influence Coefficients 
The influence coefficients turn out to be related, not surprisingly, to the 
geometry of the airfoil and the manner in which the panels are formed. Specifically, as 
derived in {Ref. 2], the A’s and B’s influence coefficients,* due to uniformly distributed 
source or vorticity are functions of: 


e The natural logarithm of the ratio a distance from - i panel control point 
(the field point) to the (j+ 1)" and j“ nodes of the j“ panel where singularities 
are distributed. 


'C’s coefficients will be needed only for unsteady flow. 
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e The angle, in radian, subtended at the i“ panel control point (the field point) by 
ene l= De and j'® nodes of the j panel where singularities are distributed. 


e The trigonometry angies of the i and j" panels. 
Referring to the geometrical quantities indicated in Figure 2.4, the 


expressions” for these influence coefficients are : 





; iar 
AY. = == ([sm(Gi oman —LiTl + cos(@.-6.)B..], forizj 
J 2X 1 J ij l J 1) 
I * ° 
= : fori = j (eqn 2.12) 
a 
AL = —{sin(@,—6.) B,. — cos(@.-@.) In 441}, fori =} 
J T 1 J 19 1 ] r.. 
ij 
= 0 ior! = j (eqn 2.15) 
| Cit — — 
BY a | oS ae L — sint Se — One |. [Or then 
J 2a l J Tx : J UJ 
i 
ma” , fori=j (eqn 2.14) 
: ] . Tey _ 
Bi. = === [ cos(@ = Gab. + sim( Gl — a iin eee ee 
J OT l ory i J ri 
i 
l ll : 
= > . fori=j) (eqn? 
where: 
Bytl a x) (xm; — X, 4 je i (ym, ~ Yip 7] 
a A [(xm, — x,)? Pr (ym, — y,)7] 
AI), a V(X FX, Zz 
yy ee Yea 
2 Actual computation uses AN = ~ B", and Bi = A‘ to reduce computing 
time. 
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Figure 2.4 Influence Coefficients due to Uniformly Distributed Singularities. 
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D. NUMERICAL SOLUTION SCHEME 
1. Rewriting the Boundary Conditions 


Using the concept of influence coefficients, the flow tangency conditions of 
Equation 2.10 can be expressed as, 


Ms 


mn 
[AT aq] + ¥ » BX + Veo sinfa-@) = 0, i=1,2,...,n (eqn 2.16) 
a 


1 


The Kutta condition of Equation 2.11, in terms of influence coefficients, looks 
as, 


n n 
a » [ A‘, q: | = >) BY = View cos(a — 8, ) 
n n 
= >| ACa qi | aay Ms B aj + Voo cos(a— ) (eqn, 2.17) 


The negative signs appearing on the left-hand-side of Equation 2.17 are a 
direct consequence of our definition of unit tangential vector. In other words, the 
tangential velocities on the lower surface panels downstream of the front stagnation 
point nave negative values. This veaturesin fact allows one to» predict the Taga 
stagnation point by interpolating the velocity distribution around the leading edge. 

2. Solving the Strengths of Source and Vorticity Distributions 

It is not difficult at this stage to see that if we collect the like terms in 

Equation 2.17 and expand Equation 2.16 for all i’s (1=1,2,....,.n), these equations 


constitute none other than a linear algebraic system of (n+1) equations as shown in 
the matrix Equation 2.18. 
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Sc ee wm avat! | ty by, | 
pious tS dont! . 2 2 
lomo nl 43 3 
ee : = (eqn 2.18) 
2 eae da 0, 
ER ee, . a Y b 
+11 +intl n+1 | 
es . Gg J 


Equation 2.18 is a set of linearly independent equations which can be easily 
solved by any standard linear system solver, one of which ts the well known method of 
Gaussian Elimination with Partial Pivoting. 

3. Computation of Velocity and Pressure Distribution 

nee the oF (j= 1,2,.....m) and y are solved, the velocities at all the panel 

control points can be evaluated. Only the tangential components exist since the norma! 


components are already set to zeroes by the flow tangency conditions. 


ee t : 
poe VY 112,50 (eqn 2.19) 
CO 
where : 
n n 
(VV) = DIAG a) + ¥D Bi + Voo cos(a—O,) , i= 1,2,..n (eqn 2.20) 


Substituting Equation 2.20 into the C ‘ equation (Equation 2.9), the pressure 
coefficients at the i panel control point is : 


eee rie oa 
(C,), ey) PN (eqme2.2 1) 
4. Computation of Forces and Moments 
The two dimensional aerodynamic coefficients of lift (Cp), drag (C,) and 


pitching moment (C mn) about the leading edge are computed by integration of the 


pressure distribution assuming constant C,, exists in each panel. The computation 1s 


et 


first done by integrating forces in the airfoil-fxed coordinate system followed by a 
rotation to the respective lift and drag directions along and perpendicular to the free 


stream (V,,) as follows : 


n 
C, = »y (C,); (Xe io x.) 
i=] 


n 
Ce ne Oey 


el 


nN 
Ca a (Co); UX 4 7%) xm; + (¥,4) 7) ym, J 


lot 


SF = = cosa + Cy sind 


Cp = Cy cosa — Gs sind 
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(eqn 2:22) 


(eqn 2zey 


(eqn 2.24) 


(eqn 2329 


(eqn 2.209 


Il. UNSTEADY FLOW PROBLEM FORMULATION 


A. OVERVIEW OF UNSTEADY FLOW MODELING 
1. Some Previews 

Having fully understood the Panel Methods formulation and solution for the 
steady flow problem, one could then venture into the interesting and complicated 
unsteady flow case. In this Chapter, we shall see how we could build the time- 
dependency into the Panel Methods solution which has been proven to be an useful 
and accurate tool for steady flow. The approach in the unsteady flow problem 
formulation will proceed, in general, in a manner similar to Chapter II. However, as we 
go along, we will pick up the highlights of the essential differences (also simularitv) 
between the two problems. Additional flow modeling of the vortex shedding process 
that greatly influences the numerical solution technique? will be discussed in details. 

2. Specific Unsteady Flow Model 

Recall that in steadv flow, the problem is considered solved as soon as the 
airfoil surface singularity distributions of source and vorticity q; i124) amon, are 
determined. These (n+ 1) unknowns are, however, time dependent in unsteady flow. 
We therefore introduce a subscript k as the time-step counter; that is, we postulate to 
solve the unsteady flow problem at successive intervals of time, starting from ty=0. At 
each time-step t, (k=1,2,....,.00), we represent the airfoil by surface singularity 
distributions consisting of source distribution (qe (j= 1,2,....,5n) and vorticity 
distribution ¥,. Again the source strengths vary from panel to panel but the vorticity 
Strength remains the same for all panels. 

The overall circulation I, at time-step t, is simply y, multiplied by the airfoil 
perimeter, £ Since the total circulation in a potential flow field must be preserved 
according to the Helmholtz’s theorem of continuity of vorticity, any changes in the 
circulation on the airfoil surface must be manifested by an equal and opposite change 
iN vorticity in the wake. We call this the vortex shedding process and postulate, as 
shown in Figure 3.1, that this shed vorticity takes place through a small straight line 
wake element attached as an additional panel to the trailing edge with uniform vorticity 
distribution (Ye We shall from now on refer to this panel as the shed vorticity panel, 
The shed vorticity panel will be established if its length A, and inclination ©,, to x- 


: Referring to the switch from a direct scheme to an iterative scheme. 
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Vortex Shedding at Time Step t, 


Helmholtz’s theorem 
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Figure 3.1 Extension of Panel Methods Representation for Unsteady Flow. 


axis of the airfoil-fixed coordinate system, satisfv the Helmhoitz’s theorem as follows, 


I 


eye) tb = Ly (eqn 3.1) 


Or AL Me = Pe 7 = OM 7 Ne) cane) 


where I, _, and y,_, are respectively the total circulation and vorticity strengths which 
are already determined at a time-step t, _, before t,. 

Let us project one time step ahead to t, ,.,, we allow the shed vorticity panel 
to be detached from the trailing edge and get convected downstream as a concentrated 
free vortex, with circulation A, (y,), or [,.,—1,, according to the resultant local 
velocity occurred at the center of the vortex core. At the same time a brand new shed 
vorticity panel is formed for the new time step and the process is repeated. Therefore 
the shed vorticity panel model provides exactly the desired communication mechanism 
to carry the solution from one time-step to another. 

We now return back to the time-step t, and immediately realise that as a 
result of this continuous vortex shedding, there has oeen a series of shedding processes 
occurred prior to t, that cummulated in a string of concentrated core vortices of 
aemcns (t,,—f,,), @,,-f,»), Fyig-2y.3)) - eee: and so on, forming the 
wake pattern behind the airfoil as shown in Figure 3.1 

The presence of the shed vorticity panel and the downstream resultant wake 
core vortices do influence the upstream flow in inviscid incompressible flow. In 
particular the shed vorticity panel itself depends on y, to determine its distributed 
vorticity (Y,,),, this in turn causes changes to (qe and y,. Moreover, the downstream 
core vortices that constitute the wake are convected under the influence of the free 
stream and the singularity distributions on the airfoil surface panels including the shed 
vorticity panel. The problem is thus seen to be coupled from this analytical 
standpoint. Putting this in simple mathematical terms, the algebriac system of 
equations (Equation 2.18), representing the flow tangency conditions and Kutta 
condition for steadv flow. are no longer linear because the coefficients as, nee. .cit- 
hand-side matrix are not constants anymore. They are function of q, and y instead. 
The presence of non-linearity is indeed what drives the solution scheme into an 


iterative type for unsteady flow. 
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3. Boundary Conditions 

We next investigate whether our unsteadv flow model 1s sufficiently 
represented, before we could proceed to search for a numerical iterative solution, by 
matching the unknowns with the available boundary conditions at time-step t,. Recall 
that we have introduced three more unknowns (Y,),. A, and ©, in addition to 
(qi), (j= 1,2,.....n) and y,. We have, however, so far only identified an extra boundary 
condition, namely the Helmholtz’s theorem (Equation 3.2) in conjunction with the flow 
tangency conditions at the n panel control points and the Kutta condition of pressure 
equilibrium at the trailing edge panels. Clearly we are in deficit of two additional 
conditions before attempting further endeavour to solve the entire system. Basu and 
Hancock [Ref. 3] suggested the following assumptions : 


e The shed vorticity panel is oriented in the direction of the local resultant 
velocity at the pane! mid point. 


e The length of the shed vorticity panel is proportional to the magnitude of the 
resultant velocity at the panel mid point and the step size of the time-step. 


Following these assumptions thus permits us to formulate two more boundary 


conditions as follows, 


eo s 
tanO, = w, (eqn 5.3) 
Ay = (yo te) LU fai | (eqn 3.4) 


where (U,,), and (V.), are the total velocity components in x and y directions of the 
airfoul-fixed coordinate system. 


The flow tangency conditions are still, 
((V").). = 0 ’  — feo ere ae (eqn 3.5) 


However, the Kutta condition must aow include the rates of change OF 
potential at the trailing edge panels (unsteady Bernoulli's equation) which can be 
related directly to the rate of change of total circulation. By using a backward finite 
difference approximation for this rate of change of total circulation, we express the 


Kutta condition as shown in Equation 3.6. 


$2 


G ol 
(V7 — (V9.7 = 21, - 0) k = 2(— 


= % 
Ot Ot 

= 2 pk Med (eqn 3.6) 
Thy 


B. RiGID BODY MOTION AND FRAME OF REFERENCE 

Consider a rigid airfoil executing a time-dependent motion, comprising linear 
translation and angular rotation about the leading edge in an inviscid incompressible 
medium. We can describe this arbitrary motion at any time instant t, as tne vector sum 
of a mean velocity — Vo, a time dependent translational velocity —[{U(t)i + V(t) j] 
and a rotational velocity —{2(t) where i & j are unit vectors in the airfoil-fixed 
coordinate system as shown in Figure 3.2. 

[f we continue, as in steady flow, to determine the flow with reference to the (x,v) 
coordinate svstem fixed on the airfoil, an observer sitting on this frame of reference 


sees an unsteady stream velocity, V made up by the vector sum of a mean 


stream’ 
velocity V.,, 2 time dependent translational velocity [U{t)i + V(t) j] and a rotational 
velocity {2(t). Therefore in this frame of reference, unlike the previous case where the 
airfoil is allowed to move only with constant linear velocity, the flow is still unsteady in 
that V.eam 1S time dependent. 


= Veo + ((U(t) i + V(t) f)] + Q(t) (yi — xj) (eqn 3.7) 


stream 


We redefine our disturbance potential to include the potential contributions @,, 


and Dey from the shed vorticity panel and the wake core vortices respectively. Thus : 

P=O, +O, +O, + OY (eqn 3.8) 
We then write the total velocity with respect to this frame of reference as, 

Poa ss on a? (eqnis5 <2} 


Notice that this total velocity is obviously NOT the absolute velocity with respect 
to an inertial coordinate system. Such an inertial coordinate system will be the one 


Where an observer only sees an on-coming flow of V,g with constant magnitude and 
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Figuré 3.2 Frame of Reference for Lasteddvel low. 


direction. We have to make this distinction clear because in our model on convection 
of core vortices, we break up the calulation into two steps; we first convect the core 
vortices using the resultant absolute velocity with respect to an inertiai coordinate 
system, followed by determining their positions with coordinates relative to the airfoil- 
fixed axes. 

The unsteady flow Bernoulli’s equation for the pressure coefficients on the airfoil 
surface must be written with respect to the airfoil-fixed coordinate system also. Giesing 


[Ref. 4] showed this to be written, in our notation, as: 





i P- Poo = ( epee \2 ss ( Pe Y ae 2 op (eqn ay 10) 
P UPVeo? Wee Ws viene eo 
where V. am? and V,...; are defined according to Equations 3.7 and 3.9. 


Equations 3.8, 3.9, and 3.10 can be correlated to their counter-parts in steady 


flow, namely 2.6, 2.8 and 2.9 respectively with V of Equation 3.7 replacing the 


stream 
Veo in Equation 2.8. 


C. TIME-DEPENDENT INFLUENCE COEFFICIENTS 
1. Definition of Time-Dependent Influence Coefficients 
The influence coefficients, Ap Aly BY and Bi involving the source and 
vorticity distributions described in Section C of Chapter II are still useful. Tnese are 
indeed time-independent coefficients since they are functions of geometrical parameters 
which are fixed in our rigid airfoil. Additional influence coefficients involving the shed 
vorticity panel and the wake core vortices must be defined. These coefficients need to 
be computed in each time step since their positions vary relative to the airfoil-fixed 
coordinate system. For that matter, as will be made clear later, those influence 
coefficients involving the shed vorticity panel must also be computed in every iteration 
within each time step for the same reasoning. 
a. More A’s and B’s Influence Coefficients 
Following the notations used previously in steady flow. we define. with the 
use of the k-subscript to denote time-dependency, additional influence coefficients 
involving uniformly distributed singularities of source and vorticity. They are the A’s 


and B’s coefficients : 


o (BY. 41), :normal velocity component induced at the i" panel control 
point by unit strength vorticity distribution on the shed vorticity panel at time 
t 
i: 
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—_ 
“ 


(Bi ne Dk : tangential velocity component induced at the i panel control 
point by unit strength vorticity distribution on the shed vorticity panei at time 
tL 

(A* aoe Me ; X-velocity component induced at eine shed vorticity panel 
point by unit strength source distribution on the j panel at time t,. 


(A” oe Ne : V-veiocilv COmponent induced at sine Shed vorticity pane! 
point by unit strength source distribution on the j panel at time ty. 


mid 
mid 


(Ee neat Dk : x-velocity component induced at a shed vorticity panel mid 


point by unit strength vorticity distribution on the j“ panel at time ty 


(BY Fea Dt y-velocitvy component induced at De shed vorticity pane 
point by unit stretieth vorticity distrioution on the j™ panel at time t, 


(A*, Oe : x-velocity component induced at — Center ror am nih 
pee by unit strength source distribution on the j“ panel at time ty 


(AY hi ye : v-velocity component induced at fe center of the hi” 
vortex by unit strength source distribution on the j panel at time t. 


(B*, Me : X-velocity component induced at the center of the h® 
vortex by unit strength vorticity distribution on the j panei at time t. 


(BY. : y-velocitvy component induced at the center 
vortex by Unit strength vorticity distrisution on the |” panel at 


mid 
core 
core 


COre 


of the nih 


Lime te 


(B “han t Dk : x-velocity component induced at the center of the ho 

vortex by unit strength vorticity distribution on the shed vorticity panel at 

te 

a een: : y-velocity component induced at the center of the h™ 
vortex by unit strength vorticity distribution on the shed vorticity panel at 
b. 


core 


core 
time 


core 
time 


New C’s Influence Coefficients 


The presence of discrete core vortices in the wake requires the definition of 


new influence coefficients involving point singularity. They are the C’s coefficients in 


our familiar notations: 


. (C “me th 
point by unit strength m 
* (¢ aa): : tangential 


doint bv unit strength m™ 
9 % . >] r 

(C peti) ae : <argl cee 
point by unit strength m™ 
(eo oT oe ; y- ore 
cea by unit strength m 


Od, : x-velocity 


vor by unit strength mm“ 


‘normal velocity component induced at the i” 


panel control 
core vortex at time ty 


velocity component induced at the i” 
COFCUVOTEEAMA TIAL, . 


panel control 


component induced at mid 


core vortex at time t 


the shed vorticity oane!l 
S 

component induced at the shed vorticity panel mid 
core vortex at time t. 
component induced at the center of the h™ core 


core vortex at time UL. 


© (C faye : y-velocity component induced at the center of the h™ core 
vortex bv unit strength m™ core vortex at time ty. 
2. Computation of Time-Dependent Influence Coefficients 
r 
(BY + De and (BY a+ De 
are computed using Equations 2.14 and 2.15 with subscript n+ 1 replacing j. Similarly, 
“4 
(AX +1) 


subsript i appropriately replaced. Also (A’, +p and (AY nj are calculated using 


are computed exactly the same way as B",, and Bi 
), and (A"yp )}, are calculated using Equation 2.12 with 0. setmie zero and 


Equation 2.13 with 9. set to zero and subsript 1 appropriately replaced. We do the 


same for ha \ Dt and (BY a using Equations 2.14 and 2.15 respectively and so 


on for all the rest of A’s ae Peemcocmictents., Ine Cs coeiiients will be computed 
with different expressions from those of A’s and B’s because they are the velocities 
induced by unit strength core vortex. It can be shown easily, from the geometry of 


Figure 3.5, that their expressions take on the following forms, 


oe) = - £088 ~ On) (eqn 3.11) 


im’k oT (r._) ). 


sinf®, — (0_.),] 


t = 
(Conk ave Cr (eqn 3.12) 


ae k 


(rd. = Vv ((xm,-x,,)? + (ym, -y_)7] 
xm, = ¥a(x, +X. 41) 
ym, = ALY; +y.44) 


= Xx coordinate of m™ core vortex at time ty 


pa 
I 


v coordinate of m“ core vortex at time t, 


<A 
{| 


Je es ) 


9. = arctan( 
ca 


Wocesn V 
(Oe = arctan(—-—™)_ 
XID, Yq 
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ee % x 2 
By the same token. (C°,_, _), and (C*,_), are computed by Equation 3.11] 
hy} y (ee ion 3. een 
while (C7 4 ee and (C*, _), are computed by Equation 3.12 if @, is set equal to zero 
and the subscript 1 appropriaielv replaced. 


D. NUNERICADSOLUTIONSCHEIteE 
. Ine flow Tangency Conditions 
The flow tangency conditions of Equation 3.5 can be written using the 
influence coefficients as follows, 


nN i 
a AM (De | + Ye » By? UV aeoegi)s © Sale 
j= j=! 
k-! 
+ (7) (Bo) + VEC) (Cy oP} = 0. ben (eqn 3.13) 
mu 


where (V__oom), is evaluated by Equation 3.7 at the i panel control point 
This equation. though it seems complex. savs nothing more than summing to 
zero all the velocity contributions due to individual singularity. Substituting (7/,,), from 


Equation 3.2, collecting like terms and rearranging the equation into, 


I P75 


| AN (A). a ty | (f/A,) (BY, a+ De a » BY 
] j= 
{ — (Mig aati ; Lf 7 (f/A,) Meal a Vk 


k-} 
—- DY uC.) ©. Bayi). ae. as (eqn 3.14) 
= 1 


2. The Iterative Solution Procedure 
Equation 3.14 is arranged in this form because we intend to solve 


Jk j= l.2......m) in termsvot 7,. Expanding Equation 3.14 feri* 1.2... results ie 
following matrix equation, 


(AMA) we, Te (eqn 3.15) 
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Figure 3.3 Influence Coefficients due to Point Singularities. 
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where [A] is an nXn matrix whose elements are known constants. { B }, and { C }, 
are n X 1 column vectors whose elements are known only if the shed vorticity panel at 
ty is established; that is, if A, and O. are known, then we can calculate ail the 
influence coefficients on the right-hand-side of Equation 3.14. We therefore make use 
of this idea to formulate our iterative solution procedure as follows : 


(1) Project the wake core vortices downstream according to the time step and the 
local resultant velocities at their respective centers with respect to an inertial 
coordinate system. 


(2) Compute the coordinates of these core vortices relative to the airfoil-fixed 
coordinate system due to its trme-dependent motion. 


(3) Start iteration cycle for current time step by initially assuming some guess 
values of A, and ©,. We can use, except for the first time-step, values 
obtained at previous time step. Compute then the influence coefficients needed 
In Equation salancigen ie: 

(4) Obtain (q.), in terms of 7, by solving Equation 3.15 as a linear system with 
two right-hand-sides by the same Gaussian elimination algorithm used in 
steady flow. 


(5) Calculate the tangential velocities at the trailing edge panels. all in terms of 
Ye: 

(6) Invoke the Kutta condition of Equation 3.6 (with some efforts in algebriac 
manipulation) to solve for y, since it is the only unknown in that equation. 


(7) Once ¥, is solved, (q,), are then known. We can then calculate the velocity 
components (U.,), and (V_), at the mid point of the shed vorticity panel. 


(8) Equation 3.3 and 3.4 hence enable us to update the values of A, and ©,. 


(9) Repeat the iteration cycle from steps (3) to (9) until converged values of A 
and ©, are obtained. Alternatively convergence can be set for (U,), and 
(Ve instead. 


(10) Compute the tangential velocities and disturbance potential at all panel 
control points in order to determine the pressure distribution which can be 
integrated to give forces and moments. 


(11) Compute the resultant velocities which occur at the centers of all the core 
vortices that will be convected down-stream. These velocities must be the 
absolute velocities with respect to an inertial coordinate svstem. 


Computation of Velocities 


3 


The iterative procedure mentioned in the previous subsection requires 
calculation of tangential velocities at the trailing edge panels and the absolute velocitv 
components (U,), and (V.),. They are computed differently due to the use of a 
different frame of reference. 
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a. Tangential Velocities on Airfoil Panels 


The tangential velocities [(V").], (i= 1,2,....,.n) at all the panel control points 
are calculated using the airfoll-fixed coordinate frame of reference as follows : 


n Tl 
(Vo. = LIAR Gt yd By 
j=l j=1 


eal Ore): E ra eee (Bi at Wk 


k-1 
i x [ (Cae @ i ov ie) ’ i 1, 2een (eqn 516) 
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b. Core Vortices Convection Velocities 


The resultant velocities at all core vortices are calculated using the inertial 
irame of reference fixed with respect to Vgo but resolving them into components in the 
directions coincide with the airfoil-fixed coordinate svstem as shown below : 
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Notice the use of V-, instead of V in Equations 3.17 and 3.18. Also 


suream 
the subscript h is just an index usaole for any core vortex. We can obtain (U.), and 
(VJ, i his replaced by n+ | in these equations. 
4. Disturbance Potential and Pressure Distribution 
a. Why We Need the Disturbance Potential 

Tne concept of disturbance potential m has been instrumental in the 
formulation of both the steady and unsteady flow problems. However, it has never 
gone beyond using it merely as a vehicle to understanding the superposition of simple 
flows. The disturbance potential need not be solved for at all in the steady flow 
problem formulation. This is because what one reallv is going after is the spatial 
derivative of this disturbance potential, 1.e. the disturbance induced velocity, from 
which the pressure distribution can be obtained. We have, in all our solutions so far, 
oeen successful in avoiding any disturbance potential calculation since the concept of 
influence coefficients allows us a direct evaluation of the velocity. Unfortunately, as 
can be seen in Equation 3.10, when we proceed further to compute the pressure 
distribution on the airfoil surface in unsteady flow, we are faced with the problem of 
evaluating the disturbance potential (®. or more precisely the rate of change of @, which 
we approximate by using a oackward finite difference expression. Therefore, the 
pressure coefficients at the i” panel control point can be rewritten, in terms of non- 


dimensional variables, as, 


D (Q.), a (Dy y 


Ck = i 7 (VA 7 C,, ae 
ke eel 


(eqn 3.12) 


eet); is the non-dimensional (by 


V oo) form of Equation 3.7 evaluated at the i panel control point. 


where (V'), is calculated by Equation 3.16 and (V 


We thus need to calculate at each time step, the disturbance potential at all 
the panel control points. Short of having to solve the Laplace’s equation by a finite 
difference scheme, we evaluate the disturbance ootential ‘p by integrating the velocity 
field in wo stages from upstream at infinity to the airtoil leading edge, then along the 
airfoil surface from the leading edge to each panel control point. Care must be taken 
here to include only the velocity contribution duc to disturbances. 

One important question arises, in this approach, as to what value of 


disturbance potential we should use at infinity before we carry out the line integral. We 
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must therefore analyse the behaviour of @ at infinity by examining the singularities that 
constitute the disturbance. They are the source and vorticity distributions on the airfoil 
surface and the core vortices in the wake. These singularities induce no velocity at 
infinity from the knowledge of simple flows. In other words, the disturbance potential 
( at infinity is independent of spatial coordinates. The next question we should ask 1s 
whether @ at infinity is time-dependent? Let us adopt the view-point that if we are at 
infinity looking at our airfoil and its associated wake, we simply see a point vortex with 
a total circulation I’) at time ty. We have already identified that I, remains constant 
by Helmholtz’s theorem. It only gets redistributed, as time progresses, over the airfoul 
surface and in the wake. Notice that the previous statement regarding what one would 
see at infinity said nothing about the source distributions. The source distributions 
though vary (or get redistributed ) as the time progresses, the total source strength 
necessarily remains zero at all time in order to enforce a closed contour representing 
the airfoil thickness. This is also the reason why the unsteady flow solution needs an 
additional model to handle the vorticity conservation since the source conservation Is 
already implicitly so for a closed contour to exist. From these discussions, we are 
certain that the disturbance potential © at infinity 1s an absolute constant (independent 
of time and spatial coordinates) whose value is fixed only by the initial condition we 
decide to start solving the unsteady problem. The actual value of @ at infinity is in fact 
immaterial so long as we Know it is constant because its value disappears conveniently 
as we subtract (@.),_, from (@.), in Equation 3.19. 
b. Computation of Disturbance Potential 

We begin by choosing an arbitrary straight line extending upstream to 
infinity from the leading edge of the airfoil along a direction parallel to Vq. For 
practical purposes, we set infinity at say ten chord lengths away from the leading edge 
since the velocities induced, at field points thereafter, by the disturbances are small 
enough to be negligible. This line is divided into z panels with element lengths near the 
leading edge comparable to the panel sizes used on the airfoil. However, the panel size 
1§ progressively increased to take advantage of the inversely decaying induced velocities 
at larger distances. We compute the tangential components of the induced velocities at 
the mid points of these panels using influence coefficients analogous to those used on 
the airfoil panels. Using subscript f to denote these panel mid-points, we can define 
influence coefficients (A’p (Aten eae (Bie) (Bien +; and (C‘.), and compute 


them using the same expressions for calculating the A’s, B’s and C’s coefficients used 
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before with cos®. replaced by (—cosa), sin8. replaced by (—sina) and subscript i 
replaced by f. With the help of these influence coefficients, the tangential velocity 


induced by disturbances at the f pane! mid point is : 
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valid for f=1,2,....,z. The disturbance potential at the airfoil leading edge is the sum of 
the products of the disturbance induced velocitv at each panel and the panel length. 
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iek = ~ LMM JU Xp2 74)? + pe Typ? | (ean 3.21) 


Similarly, for the line integral over the airfoil surface, we compute the 
tangential component of the disturbance induced velocity at the i™ panel control point 


using the following equation : 
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which 1s valid for 1=1,2,.....n. Performing the line integration by summation, the 


disturbance potential at the i nodal point on the airfoil is : 
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where Tet] denotes the panel length. 
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Finally, the disturbance potential at the i panel control point is, 
CAR amelie Pweg ne + (® ode mer) | ei (eqn 3.24) 


5. Computation of Forces and Moments 
The"Cp, C ; and C_, about the leading edge are calculated in exactly the same 
way as it is done for the steady flow problem by integrating the pressure distribution 
(See section D-4 of Chapter IT). 


EF. FLOW MODELING OF SHARP EDGE GUST FIELD 

The unsteady flow solution described so far can be extended to the study of 
airfoils penetrating a sharp edge gust by modifving the boundary conditions with the 
assumption that the gust front remains straight while passing through the airfoil. The 
same assumption has been used in both [Ref. 3] and [Ref. 4]. An additional model in 
(Ref. 3] using distribution of singularities along the gust front had successfully 
attempted to simulate the distortion of the gust front passing over the airfoil surface. 
It was shown that the pressure distributions, during the time when the gust front 
remained on the airfoil surface, were affected only at the neighbourhood of the gust 
front. The overall pressure upstream and downstream of the gust front stayed 
essentially the same. The distorted gust front model is not used in program U2DIIF. 
The use of the relatively simple yet sufficiently accurate model of a straight gust front 
affords the modifications to the unsteady flow solution to be confined only to the flow 


tangency conditions. That is to say, the expression of V in Equation 3.7 would 


stream 
include the gust velocity for panels that are already in the gust field during the 
penetration phase. Similarly, the computation of core vortex velocities using 
equations 3.i7 and 3.18 have the gust velocity added to the V., if the core vortices are 
already in the gust field. 

In an attempt to generalise the solution for cases where airfoils enter the gust 
field at an angle of attack, the convenient model used in [Ref. 3] by setting the 
computation to proceed, for the undistorted gust front simulation, so that the gust 


front always coincides with the nodal points is difficult to implement. At any one time, 
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if an airfoil enters a gust field at an angle of attack, the gust front would appear in 
oetween two nodes of a particular panel on one surface while the gust front proceeds 
from node to node on the other surface. We therefore further modify the flow 
tangency condition only on that particular panel where the gust front lies in between 
two nodes by taking the gust velocity on that panel to be proportional to the fraction 


or panel length partially submerged in the gust hteid. 
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Veo noe oN OF COMPUTER CODE U2ZDITF 


A. PROGRAM U2DIIF STRUCTURES AND CAPABILITIES 
1. Restrictions and Limitations 

The numerical formulations of both the steady and unsteady flow problems 
outlined in the previous Chapters are coded ina FORTRAN computer program called 
U2DIIF (See Appendix A for the program listings). The present solution methods 
treat the inviscid and incompressible flow as an approximation to the real flow so long 
as the viscous effect is negligible and the flow stays attached on the airfoil surface at all 
time. These restrictions are no strangers to any one who 1s familiar with anv other 
potential flow solution methods. Other than the implicit restrictions of potential flow 
solution, the method is entirely general in that the shape of airfoul is arbitrary and anv 
arbitrary continuous motion of the airfoil could be simulated using either the closed 
form (..e. explicit equations) or discrete data points to describe the time-historv of the 
translational and rotational velocities. 

The storage of the computer that carries out the calculations mav be the other 
limitation one should consider. The storage requirements grow rapidly with the number 
of panels (n) and the number of computation time steps(m). By far the prime 
contributor to this storage requirement comes from the massive amount of influence 
coefficients. The number of influence coefficients increases with the square of the 
number of panels (n*). Each time step increment adds (2n + m*) more influence 
coefficients due to the formation of shed vorticity. The current program fixes the 
maximum number of airfoil panels to 200 and the maximum allowable time steps is 
also 200. 

An additional constraint worth mentioning concerns the gust field simulation 
Whereby the current solution methodology is valid except in the use of the same 
pressure equation arising from the unsteady Bernoulli’s equation (See Equation 3.10). 
The fundamental assumption underlying the derivation of this equation is the 
irrotationality of the flow field. There is no doubt that the flow fields upstream and 
downstream of the gust front are irrotational. However, when one needs to obtain the 
pressure on the airfoil surface, an implicit integration is done across the gust front. A 


flow field inclusive of the gust front is rotational since the line integral of velocity in a 


47 


closed path does not vanish when the gust front is crossed. Failing the proper 
derivation of a new pressure equation applicabie to unsteady rotational flows, care 
must be exercised to regard the present method as an approximate solution to gust 
fields of weak strengths only. 
2. Current Structures of U2DIIF MAIN Program 

The overall program logic-ilow chart is as shown in Figure 4.1. The program 
first reads in the input data from filecode | and sets up the airfoil panel nodes and 
slopes. Immediately after that, the steady flow calculations are executed for the initial 
angle of attack a, according to the solution scheme described in Section D of 
Chapter II. The steady flow solution 1s included primarily to: 


e Provide the necessary initial parameters for the unsteady flow solution to 
proceed in time. In other words, the steady flow solution handles the V., and 
initial angle of attack a@. one decides to begin the unsteady flow calculation. 


e Allow the code to function directly as a steady flow solver as and when 
necessary without having to do the time consuming unsteady flow iterative 
solution and approach the steady flow as time approaches infinity. 


The program terminates once the. steady flow calculations are done if the 
program determines, based on the input data set bv user, no requirement for unsteady 
flow solutions. Otherwise the unsteady flow calculations will be activated bv selecting 
and computing the rigid body motions of the airfoil and the corresponding 
computation time-step size. Currently, all the time dependent motions are equation- 
generated, they are the positions and rates of the translational and rotational motions. 
Incorporated as case-runs within the program U2DIIF are the following motions : 

(1) Step change in angle of attack from any initial value. 


(2) Modified-ramp change in angle of attack about any pivot point from any 
initial value. 


(3) Harmonic translational motion at any angle of attack. 


(4) Harmonic rotational motion about any pivot point at any mean angle of 
attack. 


(5) Sharp edge gust penetration at any angle of attack. 
Snould one decide to generate the airfoil’s motion using discrete data points as a 
funcuon of time, the program could be easily modified. 
The computation time-step sizes for the harmonic translational and rotational 
motions are constant values determined by the frequencies (FREQ) and the number of 
computation per cycle (DTS). For the case of step change in angle of attack, the 


computation time-step size is progressively increasing, from a starting value (DTS), as 
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Figure 4.1 Flow Chart for U2DIIF Computer Code. 
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time increases. The case of the modified ramp change in angle of attack adopts 
initially a constant computation time-step size (DTS) during the transient rising of the 
angle of attack. Once the final angle of attack is reached, the computation time-step 
size is progressively increased also. Similarly for the case of airfoil penetrating a gust 
field, the computation time-step size (DTS) is constant during the period when the gust 
front remains on the airfoil surface but progressively increases once the entire airfoil 1s 
submerged in the gust field. These variations in computation time steps are to provide 
greater flexibility both in capturing transients and covering relatively large total time of 
computation without having to contend with the storage space requirements described: 
previously. These variations in time-step sizes described so far are associated with 
setting the input parameter TADJ to zero. If TADJ is chosen to be non-zero, all the 
case-runs would compute initially using the starting time-step sizes, based on DTS for 
non-oscillatory motions and FREQ & DTS for harmonic motions, and the program 
would prompt for an user choice of time step adjustment. If the answer is ves, the 
program would back-track the previous solution and recompute the current solution 
using an adjusted time-step size that is TADJ times the initial value (DTSY. This speciai 
time step variation feature gives the program added capability of allowing an 
interactive time step selection during the progress of unsteady flow computation. The 
ability to back-track and recompute the current solution using a different time-step size 
enhances the possibility of using program U2DIIF together with a viscous flow solver 
forming an Inviscid-Viscous-Interactive solution scheme which often requires such time 
step variations. 

The MAIN program performs the iterative solution procedures set out in 
Section D of Chapter III. The convergence check during the iterative solution is done 
through the user specified tolerance between successive iterative solutions of both 
(U.), and (V_),. The solution continues into the next time-step by selecting the time 
step size according to the particular case-run and projecting all the wake core vortices 
downstream so that their new positions relative to the airfoil at the new time step can 


Be cOmecuy determined. 


B. DESCRIPTION OF SUBROUTINES 
1. Subroutine BODY 
This subroutine is called by subroutine SETUP if the user selects an airfoil 
that is either a NACA XXXX or 230XX type. It in turn calls subroutine NACA4S to 
obtain the airfoil thickness and camber distributions and returns with the computed 


(x,y) coordinates of the panel nodal points. 
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2. Subroutine COEF 
This subroutine is called by the MAIN program in the unsteady flow 
calculations. It utilises, at each iteration cycle, the influence coefficients generated by 
subroutine INFL to calculate the coefficients of the matrix Equation 3.15 bv expanding 
Equation 3.14. These matrix coefficients are necessarily set up in this wav so that the 
source strengths could be solved in terms of the vorticity strength by subroutine 
GAUSS as a linear system with two right-hand-sides. 
3. Subroutine COFISH 
This subroutine is called by the MAIN program to set up the coefficients of 
the matrix system of Equation 2.18 for steady flow where the source strengths and 
vorticity strength are solved simultaneously by subroutine GAUSS as a linear system 
with one right-hand-side. The matrix coefficients are calculated using Equations 2.16 
and 2.17. 
4. Subroutine CORVOR 
This subroutine is called by the MAIN program at nearing the end of the 
unsteady flow calculations before starting a new time step. It computes the resultant 
convective velocities for all the wake core vortices with respect to an inertial frame of 
reference according to Equations 3.17 and 3.18 where all the appropriate influence 
coefficients are linked through common block from subroutine INFL. 
5. Subroutine FANDM 
This subroutine is used in both the steady and unsteady flow calculations. It 
is called by the MAIN program immediately after the pressure distribution over the 
airfoil panels are known so that it can perform the simple integration of pressure in the 
appropriate directions to give the aerodynamic force and moment coefficients of lift, 
drag and pitching moment about the leading edge according to Equations 2.22 through 
226; 
6. Subroutine GAUSS 
This subroutine is the standard linear svstem solver that emplovs the well 
<nown Gaussian elimination with partial pivoting und operates simultaneously on a 
user specified number of mght-hand-sides. It is cailed by the MAIN program in both 
the steady and unsteady flow calculations. In order to use GAUSS, the coefficients of 
the augmented matrix must be set up so that GAUSS will return the solutions 
replacing the corresponding columns of the augmented matrix that were iniually 
occupied by the might-hand-sides. The coefficient set-ups are done by subroutines 


COFISH and COEFF respectively for the steady and unsteady flow problems. 


Sal 


7. Subroutine INDATA 

This subroutine is called by the MMAIN program to read in the first three sets 
of data cards and returns to the MAIN program if [FLAG = 0. Otherwise it continues 
to read in the fourth data card as the NACA number corresponding to the type of 
airfoil and calculates the thickness parameters that will be used bv subroutine 
NACASS. 

8. Subroutine INFL 

This subroutine is the generator for all the influence coefficients that need to 
be stored and used by many subroutines associated with the unsteady flow calculations. 
It utilises the known relative geometrical parameters of the singularities to carry out 
computations based on Equations 2.12 through 2.15, 3.11 and 3.12. The MAIN 
program calls this subroutine in every iteration cycle of each time step so that the time- 
dependent coefficients can be updated as and when necessary. Time-independent 
coefficients are computed only once in the entire unsteady flow solutions. Those 
influence coefficients involving the wake core vortices are updated in each time step 
While those involving the shed vorticity panel are calculated as frequently as the 
number of iterations take to terminate a converged solution. It, however, does not 
compute and store those influence coefficients needed for the determination of 
disturbance potential (Equation 3.20) simply because they are used only once in each 
time step. 

9. Subroutine KUTTA 

This subroutine is called, in the unsteady flow calculations, by the MAIN 
program during every iteration cycle in each time step to invoke the Kutta condition 
for unsteady flow expressed in Equation 3.6. It calculates the tangential velocities at 
the trailing edge panels using Equation 3.16 in terms of the unknown vorticity strength 
that is manipulated and solved algebraicly. 

10. Subroutine NACA45 

This subroutine is called by subroutine BODY if the airfoil selected belongs to 
the families of NACA 4-digits airfoils or the NACA 5-digits airfoils of tvpe 730XX 
who share common thickness distributions with the 4-digits airfoils having the same 
thickness to chord ratio. The thickness and camber distribution data of these airfoils 
are calculated and returned to BODY. 
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11. Subroutine PRESS 
This subroutine is called by the MAIN program to calculate the pressure 
distribution over the airfoil panels after the iterative solution for the unsteady flow 
problem has successfully met the convergence criterion. I[t first computes the 
tangential velocities at all panel control points using Equation 3.16, then performs the 
disturbance potential evaluation at the current time step according to Equations 3.20 
through 3.24. Together with the disturbance potential data obtained from the previous 
time step. it calculates the pressure distribution using Equation 3.19. 
12. Subroutine SETUP 
This subroutine sets up the panel nodal coordinates for MAIN program by 
reading the 4 and 5 data sets of the input file if IFLAG = 1 is set. It skips the data 
reading if IFLAG = 0 and proceeds to set up the node distribution and call subroutine 
BODY to calculate the airfoil coordinates. The node distribution adopts a cosine 
formula in order to have closely packed panels towards the leading and trailing edges 
for improvements in solution accuracy. Regardless of how the nodal coordinates are 
obtained, SETUP determines the panei slopes and airfoil perimeter length. 
13. Subroutine TEWAK 
This subroutine is called by the MAIN program at every iteration cycle of 
each time step of the unsteadv flow calculations to compute the resultant velocity 
components at the mid point of the shed vorticity panel using Equation 3.17 and 3.18 
These velocity components are necessary to ensure the correct establishment of the 
shed vorticity panel length and orientation which governs the _ successful 
implementation of the iterative solution scheme for the unsteady flow problems. 
14. Subroutine VELDIS 
This subroutine returns to the MAIN program the velocities and pressure 
distributions for steady flow calculation using Equations 2.20 and 2.21. It also 
performs the evaluation of the disturbance potential at the panel control points. 
Though these disturbance potential data are not necessary for steady flow solution, 


chevy will be needed in the tirst time step of the unsteady flow pressure calcuiauon. 


C. INPUY BATA FOR PROGRAM! UZDHE 
Program U2DIIF reads tts input data from filecode 1. An example of the input 
data file is as shown in Appendix B for the case where the airfoil nodal coordinates are 


input by user. User could however let the program generate the nodal coordinates if 


the airfoil chosen happens to belong to the family of NACA 4-digits or 5-digits of tvpe 
230XX. To do this, simply change IFLAG to zero in the first item of the 3% set of 
data card and replace the nodal coordinates data in the 4 and 5™ sets of data cards 
by a single data card containing only the particular airfoil’s NACA number using 
Format (15). Figure 4.2 contains an itemised description of the sequential input 


variables. 


D. OUTPUT DATA FROM PROGRAM U2DIIF 

Appendix C contains a sample output data generated bv using the input data set 
shown in Appendix B. Due to the repetitive nature of output as the computation time 
progresses, only data at selective time steps are shown. The output data file begins with 
writing out what the program has read from the input data file followed by the 
computed nodal coordinates only if thev are program generated, otherwise proceeds to 
write the computed airfoil perimeter length. The next set of output data are the steadv 
flow solution parameters of distributed source strengths, vorticity strength. pressure 
and velocity distributions as well as the force and moment coefficients. The output data 
terminates at this point unless unsteady flow solution is required. It then prints, for 
each time step, the unsteady flow solution parameters similar to the previous output 
for steady flow with additional information pertaining to the rigid body motions and 
trailing wake vortices data. An explanation of the output variable names are listed in 


Figure 4.3. All output parameters are non-dimensional quantities. 
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Data Set #1 
eT LE 


Data Set #2 
el ee 


Dataset 45 
IFLAG 


NLOWER 
NOPRPE RK 


Data Set #4 
X(1) 


Data Set #5 
Y(I) 


Data Set #6 
ALPI 
DALP 
TCON 
FREQ 
PIVOT 


KGS IU 
ViGiaeT 


Data Set #7 
DEUER< 


Dizi 
PHASE 


Data Set #8 
TE 
DTS 
TOL 
WADI 





Format ([5) - | data card 
- Number of title cards to be used in Data Set =2. 


Format (20A4) - [TITLE data cards 
- Headings to be printed on output for case run identification. 


Format (315) - 1 data card |, 


- 0 if airfoillis NACA XXXX or 250XAX. 
- | otherwise. 


- Number of panels used on airfoil lower surface. 

- Number of panels used_on airfoil upper surtace (need not 
be the same as NLOWER). 

Format (6F 10.6) if IFLAG= 1 - variable data cards 

- x-nodal coordinates (divided by the chord length, c). A total 
of n+1 nodal points divided into 6 points per data card. 

Format (6F 10.6) - variable data cards. 

- y-nodal coordinates (divided by c) corresponding to the 

IFLAG = I. 


~ “fs 


Data Set #4 if 


Format(7F 10.6) - | data card 
- Initial angle of attack (AOA) tn deg. 


- Increment in AOA in deg for non-oscillatory motions. 
- Viaximum amplitude of AOA change in deg for rotational 
harmonic motions. 


- Non-dimensional rise time (Voot/c) of AOA for motion 
involving modified-ramp change in AOA. 


- Non-dimensional oscillation frequency (@c/V-)) for 
harmonic motions. 


- The length from the pivot point to the leading edge divided 
by c (postive aft) for rotational motions. 


- Magnitude of gust velocity (divided by V9) along Vo. 
- Magnitude of gust velocity (divided by VJ) perpendicular to Voy 


Format (3F 10.3) - 1 data card 


- Amplitude of chordwise translational oscillation divided bv c. 
(positive forward). 


- Amplitude of transverse translational oscillation divided by c. 
(positive downward). 


- Phase angle in deg between the chordwise and transverse 
translational oscillation with the latter as reference. 


Format (4F 10.3) - | data card 
- Final non-dimensional time to terminate unsteadv flow solution. 


- Starting ume-step size for non-oscillatory motions if TADJ=0.0. 4 
- Number of computation steps per cvcle tor harmonic motions. 
- Baseline time-step size for all motions if TADJ=0.0. 


- Tolerance criterion for.checking the convergence between 


successive iterations of (U,), and (V.), 


- Factor by which DTS will be adjusted. 





Figure 4.2 List of Input Variables. 
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TK - Time step t,. 

TKMtI - Time step t, |. 

ALPHA(T)  - Angle of attack at time t,. 

OMEGA(T) - Rotational velocity (positive counter clockwise) at time t,. 


U(T) - Chordwise translation velocity (postive forward) at time t,. 

wil) - Transverse translational velocity (positive downward) at ume t,. 

NITR - Iteration number. 

VXNW - Iterative solution of (U.,),. 

VY W - {terative solution of (V,),. 

Wake - Iterative solution of shed vorticity panel length A,. 

THETA, - Iterative solution of shed vorticity panel orientmtiort er 

GAMMA - iterative solution of the strength of vorticity distribution. 

J - Panel number. 

X(J) - x-coordinate of the mid point of j” panel. 

QJ) - Strength of source distribution on the j" panel. 

Cr) - Pressure coefficient at the mid point of j” pape. 

V(J) - oa pee! velocity _at the mid point of j” panel 
eferenced to the airfoul-fixed coordinate system. 

CD - Drae COcimetent. 

cL - Lift coefficient. 

CM - Pitching moment coefficient about the leading edge. 

M - Trailing wake core vortex number. 

xX(M) . = X-coordinate of the center of m core vortex. 

Y(M) - y-coordinate of the center of m™ core vortex. 

GHC - Circulation strength of the m™ core vortex. 


Figure 4.5 List of Output Variables. 
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Y. RESULTS AND DISCUSSIONS ON*CASE=RUNS 


This Chapter presents the results of numerous case-runs of U2DIIF code for the 
purpose of verifving the code. The various airfoils used in the case-runs are deliberately 
chosen to be the same as those airfoils where direct comparison of results can be made 


with either theoretical analyses and/or experimental data available in the literature. 


A. STEP CHANGE IN ANGLE-OF-ATTACK 
1. Case-Run Definitions 
Consider an airfoil initiallv at zero angle of attack to the free stream Vg that 
undergoes a step change in angle of attack at tp. The resulting flow should then 
provide the time-dependent information on the ouild-up of aerodynamic forces and 
moments on the airfoil resembling the classical results of Wagner [Ref. 5] calculated 
aased on linearised theory. Although Wagner prescribed a slightly different initial 
condition in that the aurroil is initially at rest and impulsively started at an angle of 
attack and velocity Vo, the difference is insignificant, especially for a symmetrical 
airfoil. This is because the seemingly different initial conditions when translated into 
the mathematical model means that the step change in AOA uses non Zero initial 
circulation Y, at t) with non Zero initial disturbance potential at infinity if the airfoil is 
cambered. For a symmetric airfoil, these initial values are all zeroes and therefore 
mathematically would be the same as the initial conditions prescribed by Wagner. 
2. Results and Discussions 
a. Von Mises 8.4% Thick Symmetrical Airfoil 
A 8.4% thick symmetrical Von Mises airfoil 1s used for this case-run where 
the airfoil performs a 0.1 rad (or 5.73°) step change in AOA. Figure 5.1 illustrates the 
changes in the pressure distributions over the airfoil at time instances corresponding to 
the airfoil having traveled distances. in terms of chord length. of 0.2. 0.5. 1.0. 2.0 and 
©. The associated trailing wake patterns at these time instances (less t= ©) are shown 
in Figure 5.2. The time-dependent build-up of aerodynamic coefficients of lift, drag, 
pitching moment and the circulation strength over a computation period of two 
traveled chord length are shown in Figure 5.3. Notice that the lift, pitching moment 
and circulation results are normalised by the respective steady state values at the same 


AOA. The apparently large initial loading on the airfoil shown in Figure 5.3 correlates 
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consistently with the results of the Piston Theory of [Ref. 6] which predicts the starting 


load on an arbitrary wing to be 


C; = 4a/M 
starting 
where M is the Mach number. In the case of an incompressible flow (M =0), the initial 
loading would be infinitely large. The same large initial loading was obtained by Kim 
and Mook [Ref. 7] who used continuous vorticities as panel singularities instead of our 
source and vorticity approach. Perhaps what remains most surprising is that the work 
of Basu and Hancock [Ref. 3] did not predict this initial loading, although they used 
the same singularity distributions as U2DIIF code. The initial large loading in lift and 
pitching moment decreases rapidly over a short time span, whereby the airfoil traveled 
approximately one-tenth chord length, before rising in a manner parallel to the Wagner 
Function. The drag, however, continues to decrease monotonically after the initial 
sharp fall. The circulation rises, as continuous shedding of vorticity takes place, slowly 
from the initial condition of zero to the asymptotic steady state value as time 
approaches ©. These results, disregarding the initial large loading associated with 
incompressible flow, are in close agreement with the results of [Refs. 3,4,7]. 
b. Thickness Effects on the Wagner Function 

In order to more closely correlate the results of U2DIIF code to the 
theoretical prediction of Wagner, we performed the step AOA change calculations for a 
very thin (1% thickness) NACA 4-digit symmetrical airfoil which in reality should 
represent a flat plate. The results are plotted as shown in Figure 5.4. Shown also on the 
same Figure are the results of the 8.4% thick Von Mises airfoil and a 25.5% thick 
symmetric Joukowski airfoil. The initial loading falls off less rapidly for the case of the 
simulated flat plate as compared to other thick airfoils but the subsequent rise in lift 


follows very closely the Wagner Function. 
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C, (Sectional Pressure Coefficient) 
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x/c (Chordwise Station as Fraction of Chord) 


(2) tV eaac 
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= ().2° 
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Figure 5.1 Pressure Distributions at Various Time Instances 


Resulting from a 0.1 rad Step Change in AO.\ fora 


8.4% Thick Symmetric Von Mises Airfoil. 
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Figure 5.1 
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Figure 5.2 Trailing Wake Patterns at Various Time Instances 
Resulting from a 0.1 rad Step Change in AOA fora 
§8.4% Thick Symmetric Von Mises Airfoil. 
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Figure 5.3. Time-Dependent Aerodynamic Parameters 


Resulting from a 0.1 rad Step Change in AOA for a 


8.4% Thick Symmetric Von Mises Airfoil. 
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B. MODIFIED RAMP CHANGE IN ANGLE-OF-ATTACK 
1. Case-Run Definitions 

The case of a step change in AOA can be considered as an usefui check for 
U2DIIF code since a handful of results from other theoretical analyses are available. 
However, a step change in AOA 1s practically not realisable since ail motions. short of 
having infinite velocities, take place over a finite time span. The step change in AOA 
that is practically possible is in fact some form of ramp rise over a short time span with 
large velocity. Even so, due to the inertia of the airfoil, an exact ramp rise in AOA 
does not physically describe the actual motion of the airfoil since finite time is also 
involved before the airfoil could build up its ramp velocitv. Same argument holds at the 
end of the ramp rise before the airfoil could stop at the final value of AOA. Therefore 
a so called modified ramp, with some form of rounding at the two ends of a ramp, is 
more likely to describe anything close to what 1s physically achievable. The theoretical 
work of Homentcovschi in [Ref. 8] considers the case of a flat plate that moves with 
constant velocity and changes the incidence about the mid chord, through a particular 


ramp fashion. described mathematically as, 


0° forte 
a(t) =<Sa (3 — 2t/t) t/t? for0 St St 
6a fort >t 


where 6@ is the magnitude of the AOA change and 7 is the rise time for the AOA to 
reach its final value. This particular function, plotted as shown in Figure 5.5, does in 
fact describe such a modified ramp. 
2. Results and Discussions 
a. Flat-Plate Case-Run 
Since the results of [Ref. 8] serves as another excellent source for tie 
verification of U2DIIF code, the obvious thing to do is to use U2DIIF to compute for 
the case of a flat plate. again simulated bv the 1° thick NACA-0O00L airfoil. executing 
cms modified ramp rise of 0.1 rad AOA over a rise time of 1.5 chord length. [his rise 


time is chosen simply to facilitate a direct comparison of results to {Ref. 8] which used 
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a rise time of 3 half-chord lengths. The results of computation are shown in [Figure 5.6 
and 5.7. Figure 5.6(a) takes a close look at the build up of lift during the initial period 
when the airfoil moves a distance of six chord lengths. The lft initially rises to about 
82% and then decreases to about 66% of the steady state value during the transient 
rise time. Thereafter it increases monotonically in a manner parallel to the Wagner 
Function. Figure 5.6({b) is a zoom view of the rather slow convergence of lift to the 
steady state value. It takes the airfoil to cover a distance of around 50 chord length 
before the lift builds up to almost 99% of the steady state value. The same results were 
obtained in the theoretical analysis of [Ref. 8]. Figure 5.7 shows a collection of the 
time-dependent aerodynamic parameters resulting from this particular case-run. 
b. Thickness Effects 

The same modified ramp function is used on the 8.4% thick Von Mises 
airfoil. The resulting lift-history plotted in Figure 5.8 shows a lower peak value of lift 
during the transient AOA rise as compared to the case of a flate plate though a similar 
trend of lift rise is obtained. Figures 5.9 and 5.10 show the results of pressure 
distributions and trailing wake patterns at various time instances. One could directly 
compare these Figures to the corresponding Figures arising from the step AOA change 
calculations and see the remarkable differences in transient characteristics as a result of 
varying the prescribed motions. Incidentally, one should realise that the non- 
dimensional rise time of 1.5 chord length is a deceivingly large number. In fact, when 
one converts this to the real time for an airfoil of 10 ft chord length moving at a low 
Mach number of 0.2, the rise time is indeed only of the order of 0.06 sec. which for 
practical purpose is close enough to a step. Neverthless, the transient part of the lift 


response is entirely governed by how one prescribes the transient motion. 
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Modified Ramp AOA (6a = 0.1 rad, t = 1.5) about 
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C. TRANSLATIONAL HARMONIC OSCILLATION 
1. Case-Run Definitions 
Although the U2DIIF code is capable of computing unsteady flow solution 
for any general translational harmonic motion described by a chordwise and a 


transverse components bearing a given phase relationship, 


h(t) 
h(t) 


ohy sin(@t) 
Oh, sin(t + A) 


where © is the oscillation frequency, » is the phase angle between the two oscillation 
components and 6h, & oh, are the magnitudes of chordwise and transverse oscillations 
respectively. The case-run to be presented in this section selects the motion to consist 
of onlv the transverse component, i.e. the heaving or plunging motion. A NACA-9015 
airfoil is chosen for the case-run. The airfoil is initially at zero AOA with the 
freestream Vo, and performs the plunging oscillation at an amplitude of ohy and a 
reduced frequency of @c/V 
2. Results and Discussions 

Figures 5.11 and 5.12 present the results of an airfoil executing a plunging 
motion at an amplitude of 0.018c but with two different reduced frequencies of 4.3 and 
17.0 respectively. These numbers are chosen to coincide with those numbers used in 
[Ref. 4]. Excelient correlations are obtained. Notice from these Figures that the 
oscillation frequency has a great influence on the magnitudes of the aerodynamic 
parameters due to the formation of significantly different trailing wake patterns for the 
same oscillation amplitude. Also to note is that the width of the resulting trailing wake 
is much larger than the amplitude of the oscillation, reinforcing the fact that the 
unsteady flow is strongly governed by the shed vorticity in the trailing wake. The lift 
and pitching moment oscillate at the same frequency as the airfoil motion but slightly 
out of phase, the phase differences vary with the oscillation frequency. The drag is 
nowever oscillating at about twice the frequency of the airfoil motion with a negative 
mean value. indicating that the plunging action indeed generates some propulsive 
thrust. [The same conclusion was arrived at in the experimental work of Halfman 
[Ref. 9] using a symmetrical NACA-0012 airfoil. 
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D. ROTATIONAL HARMONIC OSCILLATION 
{. Case-Run Definitions 

The case of an airfoil oscillating harmonically in pitch will be uniquely defined 
only if a pivot point is prescribed. If a fixed pivot point is used in the calculations, take 
for example the leading edge. any pitching motion about a pivot point other than the 
leading edge would need to be described as composed of a pitching and a translation 
motions about the leading edge. Program U2DIIF handles such conversion 
automatically without having the user to figure out the combined motion. This applies 
also to the case of modified ramp rise in AOA. The harmonic pitching oscillation is 


described bv, 
a(t) = 6a sin(@t) 


where 0a and ® are the amplitude and frequency of oscillation respectively. 
2, Results and Discussions 

The results for the case of the &4% thick Von Mises symmetric 2iriam 
oscillating at an amplitude of 0.01 rad (or 0.573°) at a rather high reduced frequency 
of wc’ V +5 = 20.0 about the leading edge, are shown in Figure 5.13. The lift, drag and 
pitching moment time-historv as well as the trailing wake patterns are very much 
simular to the case of a plunging airfoil at frequency of the same order of magnitude. 
The differences are clearly in the magnitudes and phase angles. These results check 
closely to those of [Ref. 3]. Figure 5.14 shows the results of the same Mises airfoil 
performing another harmonic oscillation at a lower reduced frequency of 0.8 and a 
large amplitude of 0.35973 rad about a pivot point 0.5 chord length ahead of the leading 
edge. [Ref. 4] conducted the same analysis for this pitch oscillation although the 
reason for using such a high amplitude of almost 23° was not clear. It is envisaged 
that such high amplitude may result in flow separation. Nevertheless, the case-run 1s 
carried out assuming validity of attached flow for the sake of comparing the results. 
Perhaps an innerent advantage. tn the light of U2ZDIITF code verification. with chemiiam 
or wigh amplitude in this case-run is that a discrepancy. if any, would show up 
significantly. Due to the use of different computation time-step size the pressure 
distributions on the airfoil, shown in Figure 5.14, do not correspond one-to-one at 
exactly the same angular positions as those presented in (Ref. 4]. However, the angular 
positions are matched to within 0.001°, 0.05° and 0.8° respectively for the three 


pressure distributions shown. They correlate very well to the results of (Ref. 4]. 
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Figure 5.13 Harmonic Pitching Motion about the Leading Edge 
of a 8.4% Thick Von Mises Airfoil 
6a = 0.01 rad, @c/VQy = 20.0. 
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(c) Trailing Wake Pattern atthe end of gacrcles. 
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Figure 5.14 
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E. SHARP EDGE GUST FIELD PENETRATION 
1. Case-Run Definitions 

The case of an airfoil penetrating a sharp edge gust field can ve computed 
using the U2DIIF code by specifying the components of gust velocity along and 
perpendicular to the freestream V.=, and the angle of attack of the airfoil. The results 
that are presented here consider the case of airfoil penetrating a sharp edge vertical 
gust at zero AOA, in view of generating information on the time-dependent lift 
resembling the classical results of Kussner [Ref. 10] based on linearised theorv. The 
gust front is taken to be at the leading edge at t, with only the transverse (vertical) 
component of 0.25V a. 

2. Results and Discussions 

Figures 5.15 and 5.16 demonstrate the variation of pressure distributions and 
trailing wake patterns respectively during and shortly after the gust front moves past 
the airfoil. [t 1s interesting to see that the resulting wake pattern after the entire airfoil 
is submerged in the gust field is as if being split by the gust front into two portions 
curling in Opposite directions. {Ref. 3] predicted a similar oehaviour for the case of an 
undeformed gust front. Due to the use of the modified flow tangency condition to 
handle the situation when the gust front happens to fall in between two nodal points, 
the pressure distributions, predicted by U2DIIF code, lie in between the results of the 
undeformed and deformed gust front models used in [Ref. 3]. We therefore conclude 
that this modified flow tangency condition produces sufficiently accurate results 
without adding complication in deformed gust field modeling which in turns limits the 
application to only sharp edge gusts. The present method would therefore preserve the 
generality for extending calculations to other types of gust front. Another comparison 
is made, as shown in Figure 5.17, by plotting the build-up of lift as a function of 
distance traveled by the airfoil in chord lengths. Shown in the same Figure are the 
Kussner Function and the results obtained from another case-run using a 25.5% thick 


Joukowski airfoil. 
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Figure 5.16 Trathng Wake Patterns at Various Time Instances Resulting 
from a §.4% Thick Von Mises Airfoil Penetrating a 
Vertical Sit@rp EGgevGtite of 0725 ag. . 
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Figure 5.16 . (cont’d.) 
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Figure 5.17 Normalised Time-Dependent Lift Cp. Cp due to 
Airrous of Various Thicknesses Penetrating 
a Vertical Sharp Edge Gust of 0.25V ~. 
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Vi. CONCLUDING REMARKS 


A. pre CONMENTS 

The U2DIIF computer code has been developed for the purpose of 
demonstrating the successful extension of the well known panel methods, which have 
been used extensively for steadv flow problems, into a powerful numerical tool for 
solving the unsteadv flow problems. The mathematical modeling of the various types 
of unsteady flows has been done with the goal of preserving the generality of the 
methods. The intention is to present a method that has the minimum inherent 
limitations and restrictions so that its usage for future applications and developments 
could remain appealing. 

The validation of the U2DIIF code has been done through the various case-runs 
of the numerous types of unsteady flow problems. The results of each case-run has 
been shown to be well correlated to the results obtainable from the literature in the 
form of theoretical analvses. numerical calculations based on different variants of panel 
singulariues and in cases where limited experimental data are available. The ability of 
those case-runs using an airfoil as thin as 1% thickness to produce results that 
correlate accurately to the theoretical flat plate results 1s perhaps a remarkable 


robustness possessed by the present unsteady flow solution methods. 


B. ENHANCING U2DIIF PROGRAM’S CAPABILITY 

It has been noted in Chapter IV that the current U2DIIF code limits the total 
number of panels to 200 and the total computation time steps to 200 also. These limits 
are not at all rigid and can be easily increased if the computer storage space is not 
critical. A point to note is that the computation time will grow rapidly as these limits 
are raised. The current linear system solver, the Gaussian elimination algorithm, used 
in the code must be concurrently improved upon to efficiently reduce the computation 
Pine sredmMibeaw(Or the iterations im each time step. A ciose examination of the matrix 
Equation 3.15, where the linear system solver is needed in every iteration, reveals that 
the coefficients of the left-hand-side matrix [A] are time-independent constants. 
Therefore the Gaussian elimination algorithm need only be done once for the entire 
unsteady flow calculations as far as the left-hand-side matrix 1s concerned. One could 


then perform, for each iteration, the manipulation of the two right-hand-sides 
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according to the steps taken for the reduction of the left-hand-side. This should cut 
down the computation time significantiv against the current method of manipulating 
both the left- and right-hand-sides simultaneously in each iteration. The savings in 
computation time was not pursued in the development of U2DIIF code because the 
orime concern was to demonstrate that the basic iterative solution scheme works for 
the unsteady ‘low problems. 

Another improvement that one may consider is to enable the code to be 
continuable from a time step where previous calculations were terminated. One sees 
this requirement necessary not only in the case of premature termination of 
computation due to some unforeseen circumstances, but also if one needs to prolong 
the computation time. Certainly with the current code structure, one has no choice but 
to perform the calculation nght from the beginning. 

A farther extension of U2DIIF code to the computation of unsteady flow 
problems involving multiple airfoils may be worth pursuing. Other research works that 
could be done based on U2DIIF code are in the area of incorporating more variety of 
rigid body motions into the code either in the form of closed-form equations or 
tapulated time history of the positions and rates of motions. It is important that one 
should use as close as possible, in the code, a rngid body motion that describes the 
physical motion before generating any numerical unsteady flow results for meaningful 
comparison to test data. This fact has been well illustrated and emphasized in the 
comparison of results of case-runs involving a step change to that of a ramp change in 
AOA with a fast rate which one could regard as a step in reality. However, remarkable 


difference in transient aerodynamics has been shown. 
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APPENDIX A 
U2DUF PROGRAM LISTINGS 


cecceeccececcccceccecccecececececcceccececcececeecececeeeeececeececececccc’e 


8 
C 
c 
C 
Cc 
ic 
C 


AAA 


C 
c 


aA 


PRO 


== 


GRAM UZDIIF 


C 
C 
UNSTEADY MOTION OF A TWO-DIMENSIONAL AIRFOIL C 
PN INGCONPRESomebe INVISCID FLOW C 
USING PANEL METHODS BASED ON THE HESS & SMITH C 

C 


MermmenemeeceecCcec CCCCCeCeceCcececCCCCCecccccceccccceccceececcecececceccecce 
CO 


MMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT, X(202) ,¥(202), 
COSTHE (201) , SINTHE(201),SS,NP1,NP2 


COMMON /WAK/ VYW,VXW,WAKE ,DT 

COMMON cee) VEWK VXWK 

COMMON /SING/ 0(200),GAMMA,QK(200) ,G 

COMMON /CORV/ CV(200) ,xXC(206) , $01 260). 7 TD ,CVVX (200) , CVV¥(200) 
COMMON /POT/ PHI(200) , PHIK(200 

COMMON /GUST/ UG(200) , VG(200) , GF, UGUST , VGUST 


COMMO 


DI 


= 


EXTV/ UE(200 
MENSION XXC(200),¥¥C(200) 


INPUT FROM FILE CODE 5 AND SET UP PANEL NODES AND SLOPES 


a eee 5925535 


WRITE 
1003 FORMAT 


1003) 
////,' DATA READ FROM FILE CODE 1',//) 
A 


1003 
{/ 
TA 
READ (1,501) ALPI,DALP,TCON, FREQ, PIVOT, UGUST, VGUST 
WRITE (s /501) ALPI,DALP, TCON, FREQ, PIVOT, UGUST ,, VGUST 
501 FORMAT (7F10.6) 
READ (1,501) DELHX,DELHY, PHASE 
WRITE (6,501) DELHX,DELHY , PHASE 
READ (1,501) TF,DTS,TOL,TADJ 
WRITE (6,501) TF,DTS,TOL,TADJ 
F (IFLAG .E ITE 1 


IF (IF (6,1005) 
1005 FORMAT ty COORDINATES OF AIRFOIL NODES', 


WRITE 


f 


; YC) 
TFLAG. .EQ. °35, — (6,1010) (X(I),Y(1) ,I=1,NODTOT+1) 
6,1020 


"Te 
1010 FORMAT i5,1626) 88" 


1020 FORMAT ( 


) s 
/,' AIRFOIL PERIMETER LENGTH = ',F10.6,/) 


STEADY FLOW CALCULATION AT ALPI 
ALPHA = AL 
WRITE 16 43 103 55 ALPHA 
1030 TE ta af See FLOW SOLUTION AT ALPHA = ',F10.6,/) 
IF (ALPHA .GT ‘ GO TO 200 
cOSATE = e68 (abaianer/130, 
SINALF = SIN(ALPHA*PT/120., 


CALL COFISH(SINALF,COSALF) 
CALL GAUSS(1) 

CALL VELDIS (SINALF, COSALF) 
CALL FANDM(SINALF , COSALF) 


INITIALISATION FOR UNSTEADY FLOW CALCULATION TO BEGIN 
HX = 0.0 
HY = 0.0 
HXO = 0.0 
HYO =e OFC 
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MMO 


AAN 


MO 


AAN 








DHX =0.,0 
DHY =" 000 
UX — 0) 8 
UY = O90 
ALP = ALPI 
DA =O. 0 
COSPA = 1.0 
SINDA = 0.0 
OMEGRY s 0.0 
XGF = 10) 210 
ANGLE = aE eo + ATAN(VGUST/(1.+UGUST) ) 
COSANG = 05 | ae 
SINANG = SIN(ANGLE 
DO 100 IG = 1,NODTOT 
UG(IG) = 0.0 
100 VG(IG} = 0.0 
PHA = PHASE*PI/180. 
VXW = COSALF 
VYW = SINALF 
GAMK = GAMMA 
=O 
M = 3 0 
TOLD = 0 0 
RIGID BODY MOTIONS OF AIRFOIL 
IF (FREQ .NE. 0.0) GO TO 1 
IF (DAES .BO. GeO) mcenles 
IF (TCON .NE. 0.0) GO TO 3 
ALPHA = ALPI + DALP 
COSALE = = SIN{ALPHAYPI/180. 
SINALF = SIN(ALPHA*PI/180 
3 DT = Das 
a0 = Dis 
GO TO 60 
2 ae ((UGUST . ae: 0.0) .AND. (VGUST JO. 1080)) GO TO 200 
TD = aba 
GO TO 60 
1 DI = 2.0*PI/(FREQ*DTS) 
i TD — piib 


Z = DT 
WRITE (6,1051) ; 
1051 FORMAT | LFF KKAKKKKAKKKKKAKAKKKKKKKKKKAKAKKAKKAKEKS | , / ; 
+ § KKK 


GIN UNSTEABY FLOW SORUITON = 7, 
+ | KAKKKHRRKRKKAK RAR AK RK ARK KK AKKRARARKRKKK HK I § 


40 M =M+ti1 
IF (T .GT. "F) Ge Te 200 


STORE CORE VORTEX COORDINATES FOR TIME STEP ADJUSTMENTS 
IF (M EQ. 1) GO TO 50 
M- 
51 ih =i) 
50 IF NE. ) GO me. 11 
mee (Dkr .20. 0.0) GO mo 22 
DO) GO Bo 32 
STEP CHANGE IN AOA 
IF (TADJ_ NE. C0) igomion 7 0 
TD LOAT (M+1)*DTS 
GO TO 70 
MODIFIED RAMP CHANGE IN AOA 
33. «OIF (T GE. TCON) | GO. 10 
DA = DAL 


34 
* (3.-2.*T/TCON) *(T/TCON) **2 
ALPHA = ALPI + DAL 


106 


OaAQA 


M6110 


34 


oe 


nZ0 


£21 
FLO 


iL 


COSALF = COS(ALPHA*PI/180. ) 
SINALF = SIN(ALPHA*PI/180.) 
DA = ALPHA - ALP 

COSDA = pa 189.) 
SINDA = SIN(DA*PI/130.) _ | 
OMEGA = - (DALP*PI/180.) * (6.*T/(TCON*TCON)) * (1.-T/TCON) 
DHX = pIVOT * (1.-COSDA) 
DHY = - PIVOT * SINDA 
UY = PIVOT * OMEGA 
MTCON =M 

ero 70 
DAL =). 

_ ALPHA = ALPI + DALP 

COSALF = Peete 82) 189°) 
SINALF = SIN(ALPHA*PI/180. 
DA = 16) 0 

Cesph = 1.0 

Sipe = 0.0 

OMEGA = 0.9 

DHX = One 

DHY =m). 0 

=O 


UY 
if (TADT .NE.40.0) GO TO 70 

TD = FLOAT(M+1-MTCON)*DTS 
conto. 70 


SHARP EDGE GUST (UGUST AND/OR VGUST) 


XGF = 

DO 110 IG = 1,NODTOT 

Ba = 0.0 

Pec) =0.0 | 

nC a ei 1G)=COSaLF + Y(IG)*SINALF 
XGP1 = X(IG+1)*COSALF + ¥(IG+1)*SINALF 
IF (IG .LT. NLOWER+1) GO TO 120 

ee ner .LE. XG) GO TO 110 

[einer .CE, XGP1) GO TO 111 

FAC = (XGF - XG)/(XGP1 - XG) 
pore = UGUST*FAC 

VG(IG) = VGUST*FAC 

EOuTomia/0 


Mumoscen JLE. XGP1)CO.TO 476 
Meeexcr .GE. XG) Gomroti?1 


FAC = (XGF ~ XGP1)/(XG - XGP1) 
ee = UGUST*FAC 

VGC = VGUST*FAC 

GO7TO 11 


0 

UG(IG) = UGUST 

VG(IG) = VGUST 
CONTINUE 

IF (XGF .LE. COSALF) MGUST = M 

Pee (tees .NE. 0.0) GO TO 70 
a eae COSALF) TD = FLOAT(M+1-MGUST)*DTS 
SS ~ 


TRANSLATION HARMONIC OSCILLATION 


Peer eNews 0.0) GO TO 12 
Hx DELHX eee tnt + PHA) 


HY = DELHY * SIN(FREOAT) 

DHX = HX - HXO 

DHY = HY - HYO 

UX = racemes Fa) 
UY = DELHY*FREO*COS( FREOAT) 

GO TO 70 


ROTATIONAL HARMONIC OSCILLATION 
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MAO 


V0) 


Ma 


a0 
30 


L008 
1004 


TOG 


20 


300 
400 


2001 


DAL - DALP*SIN (= SREOAT ) 
ALPYA =e Dee 

COSALF = 205 (sLPHA*PT/180. 
SINALF = SIN(A ALPHA*PI/180. 
DA = ALPHA - ALP 

COSDA = ere A*PT/130. 
SINDA = SIN(DA*PTI/180. 
OMEGA = - (DALP*PT/180.) * FREQ * COS(FREO*T) 
UM = PIVOT * OMEGA 

DHX = PIVOT * (1.-cOS@m 
Dare = - PIVOT * SINDA 


TRANSFORM CORE VORTEX COORDINATES W. R. T. NEW AIRFOIL POSITION 


TF (HM Ho. 1) GO To 80 
DO 90 I=1,M-l | 
XxC(I) = axC(T) a ae * OT 
YC(i) = Yi@ier) + Case Da 
xCO = xc(t 
YCO = YC(I 
XC(I) = XCOXCOSDA - YCO*SINDA + DHX 
7C(t = XCO*SINDA + YCO*COSDA + DHY 
CONTINUE 
WRITE (6,1001) T,DT 
SORMAT (/////,' TIME STEP TK = ',F10.6,10K,'TK - TKML = ',f10.6,/) 
WRITE 6, 1004) ALPHA, OMEGA, Ux, UY 
FORMAT (/,' ALPHA(T) = ',F10.6,5%,' OMEGA(T) = pple 
+ U(T) = RTO 6 ox, V(T) = ' F10. od L 
+ 1X,' NITR VxkW VYW WAKE THETA’ i 
CALCULATE THE TRAILING =ZDGE WAKE ZLEMENT 
WARE = E ORT (VYWAVYWHVRWAV KW) DT 


LTHEae | STAN2 (VW. VW) 

COS Tae Ne = COS (THENP | 

SENTae (NPL) S= SEyCroaetie tl 

WRITE (6,1002) NITR,VXW,VYW,WAKE, THENP1 ,GAMK 
FORMAT (15,4F10.6,E14.6 

7 haess = eae: 7 Aras ee es 
YUE = Y(NP1) + WAKE*SINTHE(NP1 

CA INFL (NITR 

CALL COEF (SINALF,COSALF,OMEGA,UX,UY) 

CALL GAUSS(2 

CALL KUTTA eh oe 


ES 
ES 


CALL TEWAK (SINALF,COSALF 

“Ou = ABS(VYW - VYWK 

Tou ABS(VXW - VAWK 

ao ((TOL1 ca TOL) .AND. (TOL2 .LT. TOL)) "6® Toma0 
= WK 


VAW = VAWK 
Li SER |G oan Sl OP 
Nee NITR 


GO * 10 
WRITE (6,1011) NITR 
FORMAT (//,' CONVERGED SOLUTION OBTAINED AFTER NITR = ',13) 
"ALL PRESS (SINALF,COSAEF , OMBGA uy Gigglia) 
I= ((UGUST .20. 0.9) .AND. (VGUST .EO. 0.0)) @O FO 300 
CALL FANDM (SINANG, COSANG) 
GO TO 400 
CALL FANDM (SINALF,COSALF) 
CONTINUE 


ADJUST TIME STEP (TADJ .NE. 0.0) IF NECESSARY 
IF (TADJ Te 0.0) GO Te 35 
a 


WRITH (5,20 
FORMAT I ' DO YOU WANT TO ADJUST TIME STEP ? © + NO, 1 © YES') 
READ ( IDT 
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3 (IDT EQ. — Go TO 95 
DT 


T TOLD : DT 

WRITE (6,1006) 

1006 FORMAT (//,' BACK-TRACK COMPUTATION AND ADJUST TIME-STEP' ,//) 
0 


GO TO 
C 
c WAKE ELEMENT LEAVES TRAILING EDGE AS A CORE<-VORTEX 
e 
5 CV(M) = yee GAMK ) 
oa = K ae o Berar Nei 
Tech = Y(NPL) + 0.5*WAKE*SINTHE(NP1 
ot = VAW 
CVVY = VYW 
WRITE ie 


1052 FORMAT (/ TRAILING VORTICES DATA', 
+ aK IM", 4X, KH), 6X, 'Y(M)!' , 6X, rCTRG! hy 


900 WRITE (8 1050)’ 1 Meme ve(1),CV(T) 
1050 FORMAT(I5,3F10. 
CALL CORVOR (STNALE,, COSALF) 


C 
c RE-INITIALISE PARAMETERS FOR NEXT TIME STEP CALCULATION 
C 
Dey 30 I = 1,NODTOT 
Q(T) = QK(T) 
PHI(I) = PHIK(I) 
30 CONTINUE 
GAMMA = GAMK 
ALP = ALPHA 
HxXO = HX 
HYO = HY 
TOLD = T 
DE = TD 
T = T + TD 
GO TO 40 
Zo) STOP 


END 
Ree epee mere oo ce gee comer c ce eel CCCCECCECCECCECEECE 
G 
SUBROUTINE BODY(Z,SIGN,X,Y) 
RETURN COORDINATES OF POINT ON THE BODY SURFACE 


C G 
c c 
ec C 
C C 
c Z = NODE-SPACING PARAMETER C 
€ X,Y = CARTESIAN COORDINATES C 
Cc SIGN = +1. FOR UPPER SURFACE G 
C ol Ons OWwE Re oURr AcE C 
c G 
C C 


Be eee ee eee eee eCeCeCe ccc Cc CeeCecececececce CCCCCCCCCCCCCCCCCCCCCCeEeCCecce 
SUBROUTINE BODY(Z,SIGN,X,Y 
COMMON /PAR/ NACA,TAU, EPSMAX, PIMAX 
ie (SIGN .LT. 0.0 Z Ae Z 
CALL NACA45 (Z, Te CAMBER BETA) 
x Cai SIGN*THICK*SIN (BETA) 
CAMBER + SIGN*THICK*COS (BETA) 


~ 
RETURN 
&ND 

OC Gree eet eree 26 eeer CCCeeec CEC CC EC See CeC Ceeeeeeeccccccccccccccccce 


SUBROUTINE COEF (SINALF,COSALF ,OMEGA,UX, UY) 


SET COEFFICIENTS “OF N”-EQUS ARISING FROM FLOW 
TANGENCY CONDITIONS AT MID POINTS OF PANELS 

SOLVING THE N-SOURCE STRENGTHS IN TERMS OF THE 
VORTICITY STRENGTH (RESULTING IN 2 RHS) 

KUTTA CONDITION IS SATISFIED SEPARATELY TO OBTAIN 
inte VORTLCI Ty STRENGTH 

THIS SOLUTION METHOD IS DESIRED FOR UNSTEADY FLOW 


) 
O 
O 
0? 


QAAANNANANANANGA 
MAANMNAMANMNOMOGHM 


109 


CeccececcccccccccccCCCCCCCCCCCCCECCCCCCECCCCECEGEC CE COC C CeCe eee ec 


AA 


AAMN 


AOAOA 


END 
cecccecccccCcceccccceGe@ccCCCCCECCCCCEEGECCECECECECECCECCCREC CCCCC CSR a 
C 


C 
C 
& 
r 
C 
€ 
C 
G 


90 


100 
140 
120 


CCCCCCCCCCCCCCCECCECCCCECECC CEE ECEECCOCECCC EEC CCEECEESCEEESECCECCECCCCCCECER 


SUBROUTINE COEF (SINALF,COSALF , OMEGA, UX, UY) 
COMMON /BOD/ IFLAG,NLOWER,NUP2ER,NODTOT, X(202) ,¥(202), 

+ COSTHE (201) 48 eUERE EEE Ss /NPi NP? 
COMMON /COF/ A(201 
COMMON /SING/ 9(200), Oe 28K (200), GAMK 
COMMON /WAK/ VYW,VXW,WAKE,D 
COMMON /CORV/ CV(200) ,XC(200) ,YC(200) ,M, TD, CCVX(200) ,cCVY(200) 
COMMON /INF1/ ne oe BEN (20%, 20m Ae een aoe sue 
COMMON /INF2/ CCN(201, 200) ,CCT(201, 200) ,CYNP1(200) ,CXNP!(200 
COMMON /GUST/ 10S (200), 'VG(200) ,XGF, UGUST, VGUST 


NEQS = NODT 

NPI = NODTOT + 1 

NP2 = NODTOT + 2 
INITIALISE COEFFICIENTS 

po 90 <T =1sNepToT 

D0 90 J = 1,NP2 


Ais) «6 00 
SET LHS MATRIX A(I,J) 


DO. 120 I = 1,NODTOT 
XMID =0.5 % ; (EIT) + K(T+1)) 
YWMID «6: = «0S * «CY(T) + VCH 

B = Om0 

Own J =e NODTOT 

A(I,J) = Gs) 

B = B + BBN(I,J) 

CONTINUE 


FILL IN TRE RIGHT HAND SIDE 


pa 2 = -B + BBN(I,NP1)*SS/WAKE 

A(I,NP2) = -BBN({I,NP1)*GAMMA*SS/WAKE 

+ + ates 3% 1 +UG¢)) )*Ccosamn vali tener 

+ - COSTHE(I)* 1.+UG(I) )*SINALF+VG(I) *COSALF+UY 

+ + OMEGA* (YMID*SINTHE(I) + XMID*COSTHE (I) ) 


ADD CORE VORTEX CONTRIBUTION 


IF (M .EQ. 1) GO To 140 

MM1 M-1 

DO 100 N = 1,MM1 

A(I,NP2) = A(I,NP2) - CCN(I,N)*CV(N) 
CONTINUE 
CONTINUE 
CONTINUE 

RETURN 


SUBROUTINE COFISH(SINALF,COSALF) 


SET COEFFICIENTS OF .LINB@AR SYSTMM - Nil souaTIONs 
NM gQUS - SLOW TANGENCY AT AID POINTS “OF PANELS 
~ EQU - XATTA CONDITION AT TRAILING EDGE PBNELS 
THIS SOLUTION “MBIHOD <5 SePECTIVE BOR creme eee, 10 
ITERATION IS REQUIRED, N-SOURCE STRENGIAS BNE 1 
VORTICITY STRENGTH ARE SOLVED SITIUETANEOUSLY 


DOAOAOIA MAI AMAA ©? 


SUBROUTINE COFISH(SINALF, COSALF ) 

, COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT, X(202) ,¥(202), 
ere cae /SINTHE (201) ,SS,NP1,NP2 

” COMMON /COF/ A(201,211),KUTTA 

COMMON /NUM/ PI,PI2INV 

KUTTA = NODTOT + 1 


110 
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M1QM 


QaOQO 


(OAM 


100 


La 


ZO 


INITIALISE COEFFICIENTS 


DO 90 Pee OS 
A(KUTTA,J) See) 


DO 120 I = 1,NODTOT 
XMID = =F): + eas 
YMID =e Yipee V( T +1.) } 
fev KUTTA) = 0.0 
-- FIND CONTRIBUTION OF J-TH PANEL 
DO 110 J = 1,NODTOT 
FLOG =o 6 
FTAN = PT 
me (5 .£Q.. 1) Couto 100 
OXI = “MID - X(J) 
DXJP = XMID - X(J+1) 
DYJ = ea ama) 
DYJP . YMED - ¥(J+1) 
FLOG = ,5FALOC( ee ee NIE | (On d Dag EE Drd)) 
STAN = AtaN2 (DYIPADKI- DXIP*OYI , DKIP*DKXJ+DYIP*DYJ) 


CTINTS = COSTHE(T)*COSTHE(J) + SINTHE(T)*SINTHE(J) 


STIMTI = SINTHE(Z)*COSTHE(3) - COSTHE(I)*SINTHE(J) 
Gay = PEt (ETANSCTINT + FLOGSTIMTJ 

3 = PIZINV*(FLOGACTIMTI - FIAN*STIMTJ 

Aer KUGbA) SeenGyRUTIA) + B 

Pec r en yim (1 .LT. NODTOT))GO TO 110 


=a Nieto Pahet TOUCHES FRAILING EDGE, 
ADD CONTRIBUTION TO KUTTA CONDITION 


a” ,J) = ACKUTTA, J) - 

A(KUTTA,KUTTA) = A(KUTTA, KUTTA) Eels) 
CONTINUS 

rILL IN KNOWN SIDES 


A(I,KUTTA+1) = SINTHE(I)*COSALF - COSTHE(1I)*SINALF 
CONTINUE 


A(KUTTA,KUTTA+1) = - tren tis c Pree eect ene 
- 


- (SINTHE(1) + SINTHE(NODTOT) )*SINALF 


RETURN 


END 
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SUBROUTINE CORVOR (SINALF,COSALF) 
GCOUPUGE THE LOCAL VELOCITIES OPSCORE VORTICES 


AANA! 
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UBROUTINE CORVOR (SINALF, COSALF) 
COMMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT,X(202) ,Y¥(202), 
+ sue /SINTHE (201) ,SS,NP1,NP2 
COMMON /SING/ 0(200),GAMMA,OK(200),GAMK 
COMMON /WAK/ VYW,VXW, WAKE , OT 
COMMON /CORV/ CV(200),2C(200) ,¥C(200) ,M,TD, CCVX(200) , CCV¥(200) 
COMMON INES / sitZ{200, Z03), at 308 201 
COMMON /INF4/ CMY(200,200) ,CMX(200, 200 
COMMON /POT/ Peay ,PHIK(200 


COMMON oo 200) ,VG(200) , XGF , UGUST, VGUST 
IF (M .EO. 1) GO TO 40 
MM1 =M- 
VELOCITY COMPONENTS OF CORE VORTICES AT CURRENT TIME STEP 
UGC = 0.0 
VGC =0)-0 


sy 


DO 10 N = eiML 


ZG . 0 (M1) *COSALE + YC(N)*SINALF 
iz (XG .GT NGF) Govwtoes 
wee = Jeucu 
vce = Veuon 
5 VY = §S*2MY(N,NP1)*(GAMMA-GAMK) /WAKE+ 

~ (1.+UGC) *SINALF+VGC*COSALF 
Vx = SS*AMY(N,NP1)*(GAMMA-GAMK) /WAKE+ 

. (1, sUGC) *COSALF -VGC =SINAL 
DO-20 Of =)1  Nebaer 
VY = V¥ + are 2) FOK(3 + eat Tre 
TX = VX - BMZ(N,J)*0OK + AMY(N,J)*GAMK 


20 CONTINUE 
aDD CORE VORTEX CONTRIBUTION 
DOws0 MG °= 1a 


Iz (MC .=0. N) GO TO 30_- 
jae cat ee 


10Q 


ae VK + CMX(N,MC)*CV(MC) 
30 CONTINUE 
Ee 
Cc CCORDINATES OF CORE VORTICES AT NEXT TIME STEP 
y= 
CCVX(N) = VK 
GGiar (HM) = VY 


10 CONTINUE 
a0 CONTINUE 
50 CONTINUE 


ccececececcece CCCCCCECECECECCCECC CEC CE CEECCECCEECESECECCCECECECCCECEC Raa 
C 
SUBROUTINE FANDM(SINAL?,COSALF) 


C 
C 
COMPUTE‘ AND PRINT OUL CDeei et G 
INTEGRATE PRESSURE DISTRIBUTION BY TRAPEZOIDAL RULE G 
G 
E 


CCCCCCCCCCCCCCCCCCCCOCCCCCCCCGECCCECCEeEcCCGGeseeeeeecececeaceccCccCcaamem 
SUBROUTINE FANDM(SINALF,COSALF) 
COMMON /BCD/ IFLAG,NLOWER,NUPPER ,NODTOT, X(202) ,¥(202), 
ze COSTHE (201), /SINTHE (201),SS,NP1,NP2 
COMMON /CPD/ CP (200 
CFX 


AMAQDAMAMNO 


F 10 
CFY = 0. 0 
CM = O80 
DO 100 I = 1,NODTOT 
XMID = Bate + es} 
Yee) = .55(Y(1) +) alee 
DX = XK ae os 5 
by mY (ey) — eae 
Cax = Cr. + Gr r)=Dx 
Cay = Ge - Crit ge. 
M = CM + CP(1)*(DKaMD 4+ DY~wmD) 


C 
100 CONTINUE 
CD 


(3 F 


CFX*COSAMR +°CFEY*SINAME 

e2 CFIACOSAL = CAN"SENELYT 
WREES (G6 ,4C0C) CD, agen 

1000 FoRMad(//, CO =! 510. CL ='@y%0 2 OT =' F105) 
ne URN 


END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC Ce ee Cee ee ee eee CCC eet 
SUBROUTINE GAUSS (NRHS) 


SOLUTION OF LINEAR ALGEBRAIC SYSTEM BY 
GAUSS ELIMINATION WITH PARTIAL PIVOTING 


COEFFICIENT MATRIA 
NUMBER OF EQUATIONS 


Ul 


ANIAAIANNIN 
ANA AA A aA A” 


A 
NEQNS 


he 


NRHS = NUMBER OF RIGHT HAND SIDES 


RIGHT-HAND SIDES AND SOLUTIONS STORED IN 
COLUMNS NEQNSt+1 THRU NEQNStNRHS OF A 


maveccecgeeccececececececececceececccccccccccccccececceccccccccccccccc 
SUBROUTINE GAUSS (NRHS 
COMMON /COF/ A(201,211), NEQNS 
NP NEONS + 1 

NTOT = NEONS + NRHS 


M0101 01010 


ANMAMAN 


ec 
C GAUSS REDUCTION 
C 
DO 150 I = 2,NEQNS 
EC 
C -- SEARCH FOR LARGEST ENTRY IN (I-1)TH COLUMN 
c ON OR BELOW MAIN DIAGONAL 
C 
IM =I-1 
IMAX = IM 
AMAX = ABS(A(IM, IM) ) 
DO 110 J = I,NEONS 
IF (AMAX .GE. "ABS (A(J, IM))) GO To 110 
IMAX = il 
AMAX = ABS(A(J,IM)) 
110 CONTINUE 
C 
Cc -- SWITCH (I-1)TH AND IMAXTH EQUATIONS 
c 
IF (IMAX .EQ. IM) GO TO 140 
Doets0 J = IM NTOT 
TEMP = Nee oS) 
ey = A(IMAX,J) 
A(IMAX,J) = TEMP 


130 CONTINUE 


ELIMINATE (I-1)TH UNKNOWN FROM 
ITH THRU (NEQNS)TH EQUATIONS 


140 po 150 J = 
R ats, TH) 7A(IM, IM) 
K = I/NTOT 
150 as.KD = A(J,K) - R*A(IM,KB) 


AAAN 


C 
C BACK SUBSTITUTION 
C 
DO 220 K = NP, 
A(NEQNS ,K) = A(NEONS K) /A(NEQNS , NEQNS) 
DO 210 L = 2,NEONS 
I = NEQNS pee 
IP =I+1 
DO 200 J = IP NEO 
200 Av eve = A(I, FR} = ad J)*A(J,K) 
210 A(I,K) = A(I,K)/A(I,15 
220 CONTINUE 
RETURN 
Bh 
ccccecc Saaeemeacerscecc cccceeceeeecccccccccecccccceeeececcccccccecceecccce 
C SUBROUTINE INDATA C 
C SET PARAMETERS OF BODY SHAPE C 
C FLOW SITUATION, AND NODE DISTRIBUTION C 
C 
C USER MUST INPUT C 
C NLOWER = NUMBER OF NODES ON LOWER SURFACE C 
C NUPPER = NUMBER OF NODES ON UPPER SURFACE C 
C PLUS DATA ON BODY AND SUBROUTINE BODY C 
C 
ceecceccecceecceeceecececcececeececececceeceececcececceccecccececeeccccc 


113 


SUBROUTINE INDATA 
DIMENSION TITLE(20) 
COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT,%X(202) ,%(202), 
f COSTHE (201) ,SINTHE(201),SS,NP1,NP2 
COMMON eae NACA, TAU, EPSMAX, PTMAX 
READ (1,501) ITITLE 
WRITE (6,501) ITITLE 
Bo 10 IT = ttre 
READ Lene) TIL 
10 WRITE (6,503) TITLE 
501 FORMAT(3I5 
502 FORMAT(20A4) 
503 FORMAT(1X,20A4) 
READ (1,501) IFLAG,NLOWER,NUPPER 
WRITE (6,501) IFLAG,NLOWER ,NUPPER 
IF (IFLAG .NE. 0) RETURN 
READ (1,501) NACA 
WRITE (6,501) NACA 


IEPS = NACA/1000 
IPTMAX = NACA/100 - LO*IEPS 
ITAU = NACA - LOOOXIEPS - LOO*IPTMAX 
EPSMAX = IEPS*0.01 
PTMAX = IPTMAX*0.1 
TAU = ITAU*0.01 
TF (IZPS .LT. 10) RETURN 
PTMAX = 0.2025 
EPSMAX = 2.6595*PTMAK**3 
RETURN 

D 


EN 
CCcceccceccccceecccCccccececcc cece Cece ce cc CC CC eee ee eee CC eee Oe CC ere ere am 
L 
C SUBROUTINE INFL (NITR) 
© 
C CALCULATE INFLUENCE COEFFICIBNTS 
C 
S 


eecceececeecccecececececcececeeecccecceeeeceeccceccccececcecececeececaam 
SUBROUTINE INFL (NITR) 
COMMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT,X(202),Y(202), 
+ COSTHE (201) ,SINTHE(201),SS,NP1,NP2 
COMMON /NUM/ PI,PIZINV 
COMMON /WAK/ VYW,VXW,WAKE,DT 
COMMON /CORV/ CV(200$,xC(200) , ¥C(200) ,M,TD,CCVx(200) ,ccVy(200) 
COMMON /INF1/ AAN(201,201),BBN(201,201 Aes Ser) eae Ce 
COMMON /INF2/ CCN(201,200),CCT(201,200) ,CYNP1(200) ,CXNP1(200 
COMMON /INF3/ AMY(200,201) ,BMY(200,201 
COMMON /INF4/ CMY(200,200),CMX(200, 200 
IF ((M .GT. 1) .OR. (NITBR .GT. 6)) Go Toweiin 


ANAMANM 


C 
E AAN(I,J) : NORMAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 
c BY UNIT STRENGTH DISTRIBUTED SOURCE ON THE J-TH PANEL 
C BBN(I,J) : NORMAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 
: BY UNIT STRENGTH DISTRIBUTED VORTEX ON THE J-TH PANEL 
DO 120 I = 1,NODTOT 
XMID = S34KCT) + ieee 
YMID = see ee) 
NO 110°oL Sar Mopser 
FLOG = Lie 
FTAN = Pi 
we (J 480. GO TO 100 
DxJ = XMID - X(J) 
DXJP = XMID - X(J+1) 
DYJ = YMID - Y(J) 
DYJP = YMID - Y(J+1) 
FLOG = ,5*ALOG( DXJP*DXIP+DYJP*DYIP ) / (DKI*DXJ+DYI*DYJ) ) 
AN = ATAN2(DYJP*DKXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 


Fr 

100 CTiid. = COMER SCORES) + SINTHE(I)*SINTHE(J) 
STIMTJ = SINTHE(I)*COSTHE(J) - COSTHE(I)*SINTHE (J) 
AAN(I,J) = PIZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 


114 


led 


MAMAMAMAAN 


140 


MNANAAMNMIAMNM 


BBN(I,J) 
CONTINUE 
CONGBHUE 
[oe 

AMID 

{MID 

BO. 130 

cLOG 


AN 
CTI. = 


BBN(I,J 
AYNP1(J) 


BYNP1(J) 


ey) 
BYNP1 (J 
CONTINUE 
DO 140 
AMID 
rp 

J 

DAT 
DXJP 
DYJ 

By IP 
FLOG 
FTAN 
CErMGe 
ST J 
AAN 2} 
Bon(l,J 
CONTINUE 


Gein o ) 


Cems) 


STIMTN 
CCN(I,N) 


nud udb nb hot iat tlt 


= PIZINV*(FLOG*CTIMTJ - FIAN*STIMTJ) 


NP1 

ees eg + a ; 
foe Ll) eS YC I+ 

= 1,NP1l 

Oe 

PL 
3 AD GO. 20.135 
XMID - K(J) 

XMID - X(J+1) 

YMID = Y(J) 


YMID - YfJ+1) 

. 5*ALOG( (DXIP*DXIP+DYIP*DYJP) / (DXI*DXJ+DYJ*DYJ) ) 
BTAN2 (DYSP*DXJ- DXJP*DYJ , DXIP*DXJ+DYJIP*DYJ) 
COSBHE (T)<COSTHE(J) + SINTHE(I)*SINTHE (J) 
SINTHE(T )*COSTHE(J) - COSTHE(1)*SINTHE (J) 


se 3 = se + FLOG*STEMT J 


PIZINV*(FLOG*CTIMTJ - FIAN*STIMTJ 


eee eee TY SE IDUCZD AT MID POINT OF WAKE ELEMENT 


(NFi=ie PANEL) BY UNiT STRENGTH DISTRIBUTED SOURCE 
ON Js 22 PANEL 


; y = VEEBOCITY INDUCED AT MID POINT OF WAKE ELEMENT 


HT a 


(NP1-TH PANEL) BY UNIT STRENGTH DISTRIBUTED VORTEX 
ON J-TH PANEL 


= MN eicoscocna. = Cee ey 
=e rn OG*COSTHE( J) + FIAN*SINTHE(J 


= 1,NODTOT 
encore + kt 
5*(Y(T) + YCI+1 
NP1 
XMID - X 
XMID - Ny J+1) 
YMiD - 
Hot) 


. 5XALOG ( DXIP*DXIP+DYIP*DYIP ) / (DKIXDKI+DYI*DYJ) ) 
ATAN2 (DYJP*DXJ-DKJP*DYJ , DXJP*DKJ+DYJP*DYJ ) 
eoenie(t) scoot ts St) See) 


ints ts 
“Kl 
Co 
i 
aonae 


S INTHE (I )*COSTHE = COstHe (1) “SINTHE( J 
Bi Ziv (branes CriIMiI Jy + cries 


PIZINV*(FLOG*CTIMTJ - FTAN*STIMTJ 


> NORMAL VEEOCITY INDUCED AT MID-POINT OF I-TH PANEL 


BY UNL stReaNnGih Noth CORE VORTEX 


: TANGENTIAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 


Lot he 


BY UNIT STRENGTH N-TH CORE VORTEX 
) i 
: 0 come. S20 


a n (E41) } 
ACY(I) + Y(I#l 
MM1 


t-4 


PS loci a 


N 
SORT (DX*DX+DY*DY) 
DX/DIST 
DY/DIST 
coams (7 *COSTHN + Sa 1) ST NTH 
SINTHE(I)*COSTHN - COSTHE(1I)*SINTHN 
= -PI2INV*CTIMTN/DIST 


nhs 


AANTIMA MOA 


MQAAAMIANHNM 


ANAAaANAYN 


Caz (aaa) 


CONTINUE 
CONTINUE 
CONTINUE 


AMID 
esD 
De 2 
DR 
DY 
DIS® 
COSTHN 
SINTHN 
Cr SrerN 
SEEN 
CEN CTN® 
ccT(i,N) 


CYNP1 (N) 


~-~ 
~ 
— 


0 


CXNP1(N) 
CxN? > N 
SENPLN 


CONTINUE 
AMY (N,J) 


; 


250 


un nnndesinnea 


= -PIZINVASPRITN/OrST 


aes 
0 


rel 
tea 


Bi 


) 


MID - ZC(N) 

SORT (DK*DK+DY*DY) 

DX/DIST 

DY/DIST 

COSTHE (1 *COSTHN + Sua) 

SINTHE(I)*COSTHN - COSTHE(I 
- 2T2INV*CTIMTN/DIST 
-PIZINV*STIMTN/DIST 


: Y - VELOCITY INDUCED AT MID POINT OF SWENE ELEMane 
(NP1-TH PANEL) BY UNIT STRENGTH N-TH CORE VORTEX 


: X - VELOCITY INDUCED AT MID POINT OF WAKE ELEMENT 
(NP1-TH PANEL) BY UNIT STRENGTH N-TH CORE VORTEX 


-PT2INV*COSTHN/DIST 
+PI2INV*SINTHN/DIST 


*SINTHN 
*SINTHN 


— 
— 
=> 
= 


: Y = VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH DISTRLBSUTED SOURGH ON DRE J-Te Patel 


BMY(N,J) : Y - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH DISTRIBUTED VORTEX ON THE J-TH PANEL 
IF (NITR .GI. 0) GO TO 530 
DO 320 N = 1,MM1 
XMID = ac 
YMID = YC(N 
pO 310 J = istepmnes 
DKJ = MTD - X00) 
DKJP = MMED @ Mi a+1) 
DYJ = YMID - Y¥(J) 
DYJP = YMID - Y(J+1) 
FLOG = . S*ALOG( ( DXJP*DXJP+DYJP*DYJP ) / (DKJ*DXJ+DYJ*DYS) ) 
FTAN = ATAN2 (DYJP*DXJ-DXJP*DYJ , DRJIP*DXJ+DYJP*DYJ) 
ee = BEF led acy oe a - oa 
BMY(N,J) = PI2ZINV*(FLOG*COSTHE(J) + FTAN*SINTHE (J 
310 CONTINUE 
320 CONTINUE 
530 CONTINUE 
DO 330 N = 1,MM1 
XMID = xC(N} 
YMID = YC(N 
i = NP1 
DXJ = XMID - X(J) 
DXIP = XMID - X(J+1) 
IYI = Hed <ev(s) 
YIP = “MID - ¥(J5+1) 
=LOG = RLOG( (DKIP*DKIP+DYIP*DYUP ) / (DKJ*DKI+DYS*DYJ) ) 
FTAN = ATAN2 (DYJP*DKJ- DXJP*™DYJ , DAIP*DKJ+DYJP*DYJ) 
4 = Pigit ( Foen ess} - L0G aes 
BMY(N,J) = PIZINV*(FLOG*COSTHE(J) + FTAN*SINTHE(J 
330 CONTINUE 
CMY(N,MC) : ¥ - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH MC-TH CORE VORTEX OTHER THAN ITSELF 
CMX(N,MC) : K - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 


STRENGTH MC-TH CORE VORTEA OTHER THAN ITSELF 


116 


IF GTR .6cr Re TURN 


Be. 420 N Selby 

XMID 7 Bei) 

YMID =n 

DO 410 MC = 1,MMl 

IF (MC .EQ. N) GO TO 410 

DX ="XMID - XC(MC 

DY = YMID - YC(MC). 

DIST = SORT (DX*DX+DY*DY) 
COSTHN = DX/DIST 

SINTHN = DY/DIST 

ae etic) = =PIZINV*COSTHN/DIST 
CMX(N,MC) = +PIZINV*SINTHN/DIST 


410 CONTINUE 
420 CONTINUE 
RETURN 


END 
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SUBROUTINE KUTTA (ALPHA,SINALF, COSALF ,OMEGA, UX, UY) 
USING KUTTA CONDITION TO DETERMINE VORTICITY 


MAMMA 


Mepmeneceeeeceeeeecceececececccececccccececececccecececececececececececece 


AAaAN 
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SUBROUTINE KUTTA (ALPHA,SINALF ,COSALF , OMEGA, UX, UY) 

COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,¥(202), 
i COSTHE (201) ,SINTHE(201),SS,NP1,NP2 

COMMON /COF/ A(201,211),NEQS 

COMMON /SING/ 0(200) ,GAMMA,OK(200) ,GAMK 

COMMON /WAK/ VYW,VXW,WAKE, DT 

COMMON /CORV/ cV(2003 ,XC(200) ,¥C(200) ,M, TD, CCVX(20 

COMMON /INF1/ a eet ace ca 3 sane: 202) BERET 

COMMON /INF2/ CCN(201,200) ,cCT(201, 200) ,CYNP1(200) ,CXNP1 

COMMON /GUST/ UG(200) ,VG(200) ,XGF , UGUST, VGUST 

DIMENSION B1(200) ,B2(200) ,AA(2) ,BB(2) 


RETRIEVE SOLUTION FROM A-MATRIX 
0 I = 1,NODTOT 
) = A(I,NP1) 

= A(I,NP2) 

FIND VT AT TRAILING EDGE PANELS 


DOe1S0 K er 


50 B 


IF (K .EQ. i) 

Baik .£0. 2) I = NODTOT 

XMID = Ono * eS + X(I+1 

YMID = 0.5 * (Y¥Y(I) + Y(Itl 

VIANG = EP I coe) 
+ + ((1.+UG(I))*SINALF+VG (I) *COSALF+UY ) *SINTHE (I 
+ + OMEGA*(YMID*COSTHE(I) - XMID*SINTHE(I) ) 

aB{K) = - AAN(I,NP1)*SS/WAKE 

BB(K = VTANG + AAN(I,NP1)*SS*GAMMA/WAKE 

Derdzo 8 = 1 N@DTOT 

BACK) 7 as + ited) - BBN(I,J)*B1(J) 

BB(K) = BB(K) - BBN(I,J)*B2(J) 


K 
120 CONTINUE 
ADD CORE VORTEX CONTRIBUTION 
IF (M .EQ. 1) GO TO 100 
MM1 =M-1 


DO 110 1,MM1 

BB(K) B(K) + CCT(I,N)*CV(N) 
110 CONTINUE 
100 CONTINUE 
130 CONTINUE 


“wa 
Ol ss 


0) ,CCVY(200) 


$605 
200 


SATISFYING KUTTA CONDITION -- SOLVE FOR VORTEX STRENGTH 


ne 








EE = AA(1)*AA(1) - AA(2)*AA(2 

FF = AA(1)*BB(1) - AA(2)*BB(2) - SS/DT 

GG = BB(1)*BB(1) - BB(2)*BB(2) + 2.*SS*GAMMA/DT 
RADI = SORT(FF*FF-EE*GG) 

GAMK = cS -FF - RADI)/EE 


CALCULA ITS SOURCE STRENGUE 


DO 160 I = 1,NODTOT 
160 QK(I) = GAMK*B1(I) + B2(I) 
RETURN 


AAO 


END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCECCECCCCCEEECEECEE CECE CEE C CCE CE CECCC Seem 
. 

SUBROUTINE NACA45(Z,THICK, CAMBER, 8ETA) 


e 

€ 

eC EVALUATE THICKNESS AND CAMBER 
e FOR NACA 4- OR 5-DIGIT AIRFOIL 
e 
Ee 


CCCCCCEECCCCCCCCCCCCCCCCCCC CECE CCCECE CEC CCEEECEC CECE CEC CCC SSC CCCs eee 
SUBROUTINE NACA45(Z,THICK, CAMBER, BETA) 
COMMON /PAR/ NACA, TAU, EPSMAX, PTMAX 
et CK = 0.0 
ie (Z .LT. leo GO 100 . 
met Ck = 5. STAUR( 5296 9RSORT (2) = 2 (wi26 + mem .55a9 
~ - ZAC 2 43-92". 2OlS) ) 
100 £E (EPSMAK smo. 9) GOmLO “250 
er beer Cin 9995 GO TO 140 
1 A EM PTMAX) GO TO 110 . 
CAMBER = EPSMAX/PTMAX/PTMAK*(2.*PTMAX - Z)*Z 
DCAMDX = 2. EPSMAX/PTMAX/PTMAX® (DTMAX 7) 


AAMNANANM 


GOLmO 126 
110 CAMBER = EPSMAX/(1.-PTMAX)**2*(1. + Z - 2.*PTMAX)*(1. - Z) 
DCAMDX = 2.*EPSMAX/(1.-PTMAX)**2*(PTMAX - Z) 
120 BETA = ATAN(DCAMDX) 
RETURN 
130 CAMBER = 0.0 
BETA = O10 
RETURN 
1g0 IE MZ .Cr. Pim) GO TO 150 
W Z/PTMAX 


CAMBER = EPSMAX*W*((W - 3.)*W + 3. - PTMAX) 


DCAMDX EPSMAX*3.%W*(1. - W)/PTMAX 
GO TO 120 
150 CAMBER = EPSMAX*(1. - Z) 
DCAMDX = - EPSMAX 
Go TO 120 


END 
CCCCCCCCOECCCCCECCCCCCCCCCCCCCEECCCCGCCCUCECC CCC CeCe ce Ce ee CeCe Cre aaa 
c 

SUBROUTINE PRESS (SINALF,COSALF, OMEGA, UX, UY) 


c 

C 

C COMPUTE UNSTEADY FEOW PRESSURE DISTRIBUTION 

C AND VELOCITY POTENTIAL AT MID-POINTS OF PANELS 
e 


CCCCCCCCCCCCCCCCEECCCCCCECCCCCECCCCCCEGEGGECC CEBBEGC CBBee 220C CCCER ea 
SUBROUT Tae PRESS (3 ee eee JOMEGR Jem, Ua; 


MAAANANAN 


COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT,<(202),7(202), 
+ Soe 20). /SINTHE(201),SS,NeP1,NP2 
COMMON /CED/ CP(200 


COMMON /NUM/ PI, 

COMMON /SING/ 9 (500). "CRRA ,QK(200) , ,GAMK 

COMMON /WAK/ VYW,VXW,WA 

COMMON /CORV/ cv(200}), 80 (200). ¥C(200) ,M, TD, CCVX(200) , CCV¥(200) 
COMMON /INF1/ AAN(201,201 BEN (201, ae} eee (20) ARIE cla 
COMMON /INF2/ CCN(201, 200 CCT(201,200) , CYNP1(200) , CXNP1(200 
COMMON /POT/ PHI(200 200 

COMMON /GUST/ UG(200 "YG(200), XGF , UGUST , VGUST 

COMMON /EXTV/ UE(200 


118 


MMO 


AAAANNA QAAaANY 


AAD 


WRITE (6y 1000») 
PIND TANGENTIAL VELOCITY VT AT MID-POINT OF [-TH PANEL 


DO=730°. i = L, NODTOT 

XMID ate i ee if ~ et 

YMID = Q. + Y(I+1 

DX = Cif re1) - a 

DY = (V(T+1) - YT 

DIST = SORT (DX*DX+DY*DY) | | 

VSX = Hi OR ee SE + OMEGA*YMID + UX 

VSY = (1.+UG T) ASINALF+VG{I)*COSALF - OMEGA*XMID + UY 

VS. = VSK*VSK + VSY*VSY 

VTANG = (fi. eet Ta CT ACOSER PSUS eTNIRE CE 
~ + ((1.+UG(I *SINALF+VG(T *COSALF+UY ) *SINTHE (I 
+ + OMEGA*(YMID*COSTHE(I) - XMID*SINTHE(TI)) 

VTFREE = VTANG 

VTANG = VTANG + SS*(GAMMA-GAMK ) *AAN(I,NP1)/WAKE 

pemeee 0 = 1 ,NODTOT 

VTANG = VTANG - BEN (I, J)*OK(J) + AAN(I,7)*GAMK 


120 CONTINUE 
ADD CORE VORTEX CONTRIBUTION 
ne SO mcONTO 150 


MM1 =M-i 
DO 140 N = 1,MM1 
VIANG = VTANG + CCT(I,N)*CV(N) 


14Q CONTINUE 
150 CONTINUE 


PHIK(I) = (VTANG-VTFREE )*DIST 
aia = VS - VTANG*VTANG 
2 1 = VTANG 


130 CONTINUE 
COMPUTE DISTURBANCE POTENTIAL BY LINE INTEGRAL OF VELOCITY FIELD 
INTEGRATION FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 


NPHI = 10 * NLOWER 

PINK = 10-0 

XL = 0.0 

homeo) = L = 1, 

FRACT = FLOKT(L) /FLOAT (NPE) 

XLP = -10.0 * (1.0 - COS(0.5*PI*FRACT) ) 
DELX = xL - NLP 

XMID = Bee) Sets 

YMID = 0.5%*(XL+XLP)*SINALF 


AL XLP 
VELXK = UGUST 
ADD CONTRIBUTION OF J-TH PANEL 


Bomzor J = NPL 

DXJ = XMID - X(J) 

DXJP = XMID - xX(J+1) 

OYJ = YMID - Y(J) 

DYJP = YMID - 7(J+1) 

FLOG = =. 5*ALOG( ( DXJP*DKIP*DYIP*DY JP ) / (DXI*DXJ+DYI*DYS) ) 
FTAN = ATAN2(DYJP*DKJ-DXUP*DYJ , DXIP*DXI+DYIP*DYJ) 
CAMITI = ~COSALF <COSTHE (J) - SINALESS INTHE (J) 
SALMTJ = -SINALF*COSTHE(J) + COSALF*SINTHE (J) 

APY = PIZINVS (FTANTCALHTJ + ELOGSSALNTS) 

BPY = PI2INV*(FLOG*CALMTJ - FTAN*SALMTJ 

IF (J .EQ. NP1) GO TO 40 

VELK | = VELX - BPY*OK(J) +GAMK*APY 

Eonro 29 


40 VELX = VELX + SS*APY*(GAMMA-GAMK) /WAKE 
20 CONTINUE 


ee 


ADD CORE VORTEX CONTRIBUTION 


610 


IZ (M (=O. 1) Gommonmso 
Med = Heer f 
DO 60 N = 1,MM1 
Dx =a a 
DY = YMID - YC(N) ~ 
DIST = SORT (DX*DX+DY*DY) 
COSTHN = DX/DIST 
SINTHN = DY/DIST. | 
SALMIN = -SINALF*COSTHN + COSALF*SINTHN 
CPT = ~DI2INV*SALMTN/DIST 
60  VELX = VELK + CPT*CV(N) 
50 CONTINUE | 
PINK = PINK + VELX * DELX 
30 CONTINUE 
e 
e COMPUTE DISTURBANCE POTENTIAL AT MID-POINT OF I-TH PANEL 
C 
€ LOWER SURFACE 
C 
DO 230 I = 1,NLOWER 
PH = -PINK 
D0 240 J = I,NLOWER 
240 PH = PH - PHIX(<) 
PUIK(I) = 2H 
230 CONTINUE 
DO 270 I /NLOWER-1 
PHIK(I) = 


= 1 
0.5*(PHIK(I) + PHIK(I+1)) 
270 CONTINUE | 
PHIK(NLOWER) = 0.5*(PHIK(NLOWER) + PINK) 
UPBER SURBERCE 


I = NODTOT,NLOWERF1,-1 
= -PINKA 
J 


AAO] 


DO 250 
a 


P 
DO 260 = NLOWER+1,I 
260 PH = PH + PHIK(J 
PHIK(I) = PH 
250 CONTINUE 
90 280 I = NODTOT,NLOWER+2,-1 
PHIK(I) = 0.5%(PHIK(I) + PHIK(I-1)) 
280 CONTINUE 
PHIK(NLOWER+1) = 0.5*(PHIK(NLOWER+1) + PINK) 


COMPUTE CPeAT MGD POPNT OF I=TH PANEL 


DO 290 I = 1,NODTOT 
XMID = .5*(X(I) + X(t+1)) 
CP(I) = cPp(I) - 2 oar ee 
WRITE (6,1050) I,XMID,QK(I),GAMK,CP(I) ,UE(I) 
290 CONTINUE 
1000 FORMAT(/,4X,'J' 4X, 'X(J)',6X,'Q(J)',5X, 'GAMMA' , 5X, 
+ 'CP(J)', 6X, 'v(g)",/) 
LQ50 PORMAT(I5,6F10.6) 
RETURN 
coND 
SCCICCCCCCCCCCCECCCCCCCCCC CEG EGGCC EG CCS CCGG GG ee 6 a 6 Oe eee 


SUBROUTLAE SiemeOP <£ : 


C C 
C C 
C SETUP COORDINATES OF PANEL NODES AND SLOPES OF PANELS C 
c COORDINATES ARE READ FROM INPUT DATA FILE UNLESS ¢ 
C 
C C 
C C 


AAD 


THE AIRFOIL IS OF NACA KXXXX OR NACA 230KX TYPE 


CCCCCCCCCCCCCECCCECCCCCCCCCCCCCCCECCCCCC ESC ee COC CC eee ee ce CC eC Cee 
SUBROUTINE SETUP 
COMMON /BOD/ IFLAG,NLOWER ,NUPPER,NODTOT,X(202),Y(202), 
7 COSTHE (201) SINTHE (201) ,S5S,NPL,1eZ 
COMMON /NUM/ PI,PIZINV 


120 


PI = 3.1415926585 
PIZINV = .5/PI 
C 
C SET COORDINATES OF NODES ON BODY SURFACE 
C 
IF (IFLAG .NE. Racor TO 10 
NPOINT NLOWER 
SIGN een 
NSTART = 0 
DO 110 NSURF = 1,2 
DO 100 N = 1,NPOIN 
FRACT = FLOAT(N-1)/FLOAT(NPOINT) 
Z = .5%(1. - COS(PI*FRACT)) 
i = NSTART + N 
CALL BODY(Z,SIGN,X(I) ,Y(I)) 
100 CONTINUE 
NPOINT = NUPPER 
SIGN =e 
NSTART = NLOWER 


110 CONTINUE 


NODTOT = NLOWER + NUPPER 
ett) =. (1) 
Y(NODTOT+1) = Y¥(1 
GO TO 20 
10 NODTOT = NLOWER + NUPPER 
READ (1,501) (X(I),I=1,NODTOT+1) 
WRITE (6,501) (X(1I),I=1,NODTOT+1) 
READ (1,501) (¥(I),I=1,NODTOT+1) 
WRITE (8 501) (¥(1),I=1,NODTOT+1) 
501 FORMAT (6F10.6) 
20 NPI = NODTOT + 1 
NP2 = NODTOT + 2 
C 
C SET SLOPES OF PANELS AND CALCULATE AIRFOIL PERIMETE 
C 
SS = 
DO 200 I = 1,NODTOT 
DX = X . = z A 
DY = Y(T+l 
DIST : SORT (Dx*DK +DYADY) 
es 2 = ms 
Same (1) = DY/DIST 
HE(I) = DX/DIST 


GOST 
200 CONTINUE 
RETURN 
END 
mame oeceewmeeecO CCCCCECCCCCECCOCCCECCCCECCECC etc eteeececcceccccccccceccc 
SUBROUTINE TEWAK (SINALF,COSALF ) 


COMPUTE WAKE ELEMENT AT THE TRAILING EDGE 


mpeemeedeeeccecceccccecceccccccccecccccecccceccceeeccaeeetcccecceccccccc 
SUBROUTINE TEWAK (SINALF,COSALF 
COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT ,X(202) ,¥(202), 
e COSTHE (201) ,SINTHE(201),SS,NP1,N>P2 
COMMON /COF/ A(201,211) ,NEOS 
COMMON /SING/ 0(200) ,GAMMA,OK(200) ,GaMK 
COMMON /WAK/ UUW, VEW WAKE DT 
COMMON /WAK2/ VYWK,VXWK 
COMMON /CORV/ cV(200), RC (200) , ¥C(200) ,M, TD, CCVX(200) , CCVY(200) 


AAAMA 
ANAAN 


«16 
XMID*COSALF + YMID*SINALF 
Ls AGE GO TOMO 


COMMON /INF1/ euacot’ ,20 54 BEN (201, 201), sie 02) ENP (201) 
COMMON /INF2/ CCN(201,200) ,CCT(201,200) ,CYNP1(200) ,CXNP1(200 
COMMON /GUST/ UG(200),VG(200) ,XGF, UGUST, VGUST 
XMID = 0.5 * (X(NP1) + oes 
wad = 0.5 * (Y(NP1) + Y(NP2 
UGW = 0.0 
VGW = 0 
G 


XG 
IF (XG 


121 


UGW 
VGW 
10 VYWK 
VRWK 
DO 22z0 
VYWK 
120 VXWK 


aaAm 


DO 130 
VYWK 
130 VXWK 


= UGUST 


= (1.+UGW)*SINALF+VGW*COSALF 

= (1.+UGW)*COSALF-VGW*SINALF 

J =sep nor | 

= VYWK + AYNP1(J)*QK(J) + BYNP1(J)*GAMK 
= VXWK - BYNP1(J)*QK(J) + AYNP1(J)*GAMK 


ADD CORE VORTEX CONTRIBUTIOCEsI 
LEeGhie. 
MM1 


EQ. 1) Gomrouie 

= Me 1 

N = 1,MM1 

= VYWK + CYNP1(N)*CV(N) 
= VXWK + CXNP1(N)*CV(N) 


140 CONTINUE 


RETURN 


END 
CCccCcCccCCCCCCCCCCCCCECECCCCCCCCCECECCCE EEE CE CEC CECCC CCR e eee eee 
. 

SUBROUTINE VELDIS(SINALF,COSALF) 


S 
C 
e 
ie 
C 
C 


COMPUTE STEADY FLOW PRESSURE DIS@Ris0Tig)) 
AND VEEOCITY POTENTIAL AT HiD-POINTS OF Seaiees 


CccccccccccecccCCCECCCCCCCCECCCCCCCC EG Gee eee Ce eee eee CC CGS eee 


SUBROUTINE VELDIS(SINALF,COSALF 


COMMON /BOD/ IFLAG,NLOWER ,NUPPER,NODTOT,X(202) ,Y¥(202), 
+ COSTHE(201) ,SINTHE(201),SS,NP1,NP2 
COMMON /COF/ A(201,211),KUTTA 
COMMON /CDD/ CP(200) 
COMMON /NUM/ PI,PI2INV 
COMMON /SING/ 2600) GAMMA ,QK(200) , GAMK 
COMMON /POT/ PHI(200) , PHIK(200) 
COMMON /GUST/ UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON /EXTV/ UE(200 
WRITE (6,1000) 
C 
C RETRIEVE SOLUTION FROM A-MATRIX 
c 
DO 50 I = sep ToT 
50 Q(I) = A(I,KUTTA+1) 
GAMMA = A(KUTTA,KUTTAt1) 
¢ 
: FIND VT AND CP AT MID-POINT OF I-TH PANEL 
DO 130 & *1,NODTOT 
XMID = oF ei + X(I+1 
YMID = 5° (a1). toe) 
VTANG = COSALF*COSTHE(I) + SINALF*SINTHE(I) 
; VTFREE = VTANG 
Gc ADD CONTRIBUTION OF J-TH PANEL 
E 
DO 120 J -®i1,NODTOT 
FLOG w Onc 
STA = >T 
iP to .BO. 3) Ge @on100 
DK = KHID - K(J 
DXJP = XMID - X(J+1) 
DYJ = YMID - Y(J) 
DYJP = YMID - Y(J+#1)_ 
FLOG = .5*ALOG( DXJP*DKJIP+DYJP*DYJIP ) / (DKIXDKJ+DYI*DYJ) ) 
FTAN = Aree acostia ay" k aT DXJ+DYJP*DYJ) 
100 CTIMTJ = COSTHE(I)*COSTHE(J) + SINTHE(1I)*SINTHE(J) 
STIMTJ = SINTHE(I)*COSTHE(J) - COSTHE(I)*SINTHE(J) 
AA = PI2ZINV hear + Pe ewot nas: 
B = PIZINV*(FLOG*CTIMTJ - FTAN*STIMTJ 
VTANG = VTANG - B*Q(J) +GAMMA*AA 


120 CONTINUE 


[22 


MAAMMNAMA 


() 


MOQ 


20 


AAANDN 


AAG 


20 
30 


ANAAN 


240 
230 


AAD 


260 
250 


1.0 - VTANG*VTANG 
VTANG 
(6,2050) =, XMibeO( 1), GCAMMeeeP (I) , VE (I) 
INITIAL SeT-UP FOR DISTURBANCE POTENTIAL CALCULATION 


2) 

nj 
t4i-~ 
tint 
trj~_- 
i il 


DX = 4(T+1) - X(T) 

De = Yoer) « Y(L 

Dist = SORT (DX*DX+DYADY) 
PHI(I) = (VTANG-VTFREE)*DIST 
CONTINUE 


COMPUTE DISTURBANCE POTENTIAL BY LINE INTEGRAL OF VELOCITY FIELD 
INTEGRATION FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 


NPHI = 10 * NLOWER 

>IN = 0.0 

XL = 0.0 

ee Oe Ae = 1 

eer = FLOAT(E) /FLOAT (NPHT) 

XLP = -10.0 * (1.0 - COS(0.S*PI*FRACT)) 

DELX = XL ‘- ° eEP 

XMID = ie ee =) SCOSALA 

YMID = 0.5%*(XL+XLP)*SINALF 

ae =e ee 

VELK = UGUST 

ADD CONTRIBUTION OF J-TH PANEL 

DO 20 J = 1,NODTOT : 
DX = XMID - X(J) 

DXIP = XMID - X(3+i) 

DY = YMID - ¥(J) 

DYIP = YMID - ¥(J+1) 

FLOG = .S*ALOG( ( DXIP*DKIP+DYIP*DYIP ) / (DXI*DXJ+DYJ*DYI) ) 
FTAN = ATAN2(DYJPADXJ-DXIP*DYJ , DXJPADKI+DYIP*DYJ) 
emirsl = eo - eee gee |) 
SALMTJ = -SINALF*COSTHE(J) + COSALF*SINTHE (J 

APY = ee cree iT + ae 

BPY = PIZINV*(FLOG*CALMTJ - FTAN*SALMTJ 

VELX = VELKX - BPY*Q(J) +GAMMA*APY 
CONTINUE 

PIN = PIN + VELX * DELX 
CONTINUE 


SoUPUIE DISTURBANCE POTENTIAL AT MID-POINT OF I-TH PANEL 
LOWER SURFACE 


DO 230 I = 1,NLOWER 

PH = -PIN 

DO 240 J = I,NLOWER 
PH = PH - PHI(J) 

PHI(I) = PH 

CONTINUE 

DO 270 I = 1,NLOWER-1 

Eakin = 0. 5*(PHI(T) + PHI(I+1)) 
CONTINUE 


PHI(NLOWER) = 0.5*(PHI(NLOWER) + PIN) 
UPPER SURFACE 


bo 250 I= NODTOT , NLOWER+1 , -1 
H = - 
DO 260 J= NEOWER+1, I 
PH = PH + PHI(J) 
PHI(I) = PH 
CONTINUE 


DO 280 I = NODTOT,NLOWER+2,-1 


Ie 


t— in) 
O3 


> 


>. 


Piel (1) 


O25 0PHI (1) + Peer) 


CONTINUE 

SHi (NLOWER+1) |= ea T(NLOWER+1) + 
FORMAT (/, ax "X( 8)! 6%, Clooney 
+ ees) ' aa Jy) 
TORMAT(I5,5F10.6) 

RETURN 

=ND 
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AN 


EXAMPLE INPUT DATA FOR PROGRAM U2DiIF 


APPENDIX B 


THIS IS AN EXAMPLE OUTPUT DATA OBTAINABLE FROM PROGRAM UZDIir 


AIRFOIL 


PANEL NUMBER 


: MISES 8.4% THICKNESS (COORDINATES ARE INPUT BY USER) 


AIRFOIL MOTION 


INITIAL AOA 
FINAL AOA : 7. 
AOA RISE TIME 
COMPSGATION TIMESTEP 


: 0.05 DURING TRANSIENT MOTION, 
GRESSIVELY AFTER FINAL AOA 


: NLOWER = 25 , NUPPER = 25 
: MODIFIED RAMP AOA CHANGE ABOUT MID CHORD 
SZ oe DEGREES 
> DEGREES 
*: 1.9 CHORD LENGTHS 


INCREASES 
Po. REACHED 


es 
KAKKARKKKK] RRAKKKKKK YD RRAKAKKRKZAAKKRKKKKRAY RA KKRRAKKG RRRKKRRRAKG RRR AKRRRRT 


Ol 


° 
s 
e 
e 
° 
e 
e 
° 
td 
e 
e 
e 
@ 
° 
° 
° 
Ld 
e 


1 
Q 
0 
Q 
0 
0 
0 
0 
0 
0 
-0 
=) 
=0 
=i) 
0 
0 
0 
0 
0 
2 


25 
000000 
851308 
514918 
£73061 
003767 
W913 93 
By 208Z 
HoL7 50 
980866 
QO00000 
O17837 
W597 12 
036360 
006259 
928555 
042083 
026618 
002784 


2.90000 


000000 
000000 


DO NODDDOVVVDDD0OO0D0O000 


Zao 

.994858 
OOo 
-453098 
ee Odes 
-000000 
EZ IOLS 
S095 
soO Solo 
-994858 
.000732 
022205 
041314 
032820 
.000000 
-032820 
-041314 
-UEEZoo 
-000782 
000000 
-000000 
-050000 


© 


- 980866 


Oey os) D5 


a) 
0 


- 392084 
MOSS 
-003767 
mt SOoe 
so2915 
coos0S 
000000 
.002784 


sWiecola 
0.042083 
J0Z8555 
~006259 


-036360 


I i ' 
OODOQVOD000O DOODO0O00°O 


-958884 
-695948 
. 332794 
sWo 7146 
-015008 
TALLOoS 
paola 
893455 


Oss 21 
Wear. 
.041979 
s0Z5G01 
,OLZ2379 
039025 
-037341 
-013459 


0.0 


-900000 
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area lelt@?1°1?1 Sloe lelelelelelele) 


-eaeoso 
neo T 2) i 
27 oni 
~O3Soe8 
jUSsaou 
Bay OLD 
SS MPASTS 
,929050 


pU09S 52 
-034289 
.040979 
Ore 2218) 
sOLS220 
-040979 
-034289 
IOs o 1 


0.5 


OOODOOD9ND00 0 


=) 
-0 
= 


\ 
OOO0O 


sO952 30 
obo 
S2aa005 
.015010 
elsye gas, 
moo2 lt oe 
9246 
-958864 


~O13459 
.037341 
sJo20 06 
Oo 79 
Zao D1 
-041979 
le Won 1 
On Z 1 


0.0 


0 


.00 


APPEN DIERSC 
EXAMPLE OUTPUT DATA FROM PROGRAM U2DIiF 


DATA READ FROM FILE CODE 1 


eo 0. 0 eee 


ee eee eee ees a 


as io een 6 eee, 


THIS IS AN EXAMPLE OUTPUT DATA OBTAINABLE FROM PROGRAM U2DIiF 


AIRFOIL 


PANEL NUMBER 


AIRFOIL MOTION 
INITIAL AOA 


FINAL AOA : 7 
AOA RISE TIME 


COMPUTATION TIME STEP 


: 2.9 DEGREES 
.59 DEGREES 


, NUPPER = 25 
: MODIFIED RAMP AOA CHANGE ABOUT MID CHORD 


> 1.5 GHORD LENGTRs 


‘ . . 


: 0.05 DURING TRANSIENT MOTION, 


Loe 2a 
1.000000 0.994858 0.980866 0.958884 0.929536 0.893455 
0.851308 0.803815 0.751753 0.695948 0.637271 0.576620 
9.514918 0.453098 0.392084 0.332794 0.276105 0.222865 
0.173861 0.129819 0.091393 0.059146 0.033560 0.015010 
0.003767 0.000000 0.003767 0.015008 0.033560 0.059146 
0.091393 0.129819 0.173861 0.222865 0.276105 0.332791 
0.392082 0.453095 0.514915 0.576617 0.637266 0.695946 
0.751750 0.803815 0.851308 0.893455 0.929536 0.958884 
0.980866 0.994858 1.000000 
0.000000 -0.000782 -0.002784 -0.005721 -0.009351 -0.013459 
-0.017837 -0.022285 -0.026618 -0.030671 -0.034289 -0.037341 
-0.039712 -0.041314 -0.042083 -0.041979 -0.040979 -0.039096 
-0.036360 -0.032820 -0.028555 -0.023651 -0.018220 -0.012379 
-0.006259 0.000000 0.006259 0.012379 0.018220 0.023651 
0.028555 0.032820 0.036360 0.039096 0.040979 0.041979 
0.042083 0.041314 0.039712 0.037341 0.034289 0.030671 
0.026618 0.022285 0.017837 0.013459 0.009351 0.005721 
0.002784 0.000782 0.000000 
2.500000 $.000000 1.500000 0.000000 0.500000 0.000000 
0.000000 0.000000 0.000000 
2.000000 0.050000 0.000100 0.000000 
AIRFOIL PERIMETER LENGTH = 2.018599 
STEADY FLOW SOLUTION AT ALPHA = 2.500000 
a Q(J) GAMMA cP(J) 7(5) 
1 0.997429 0.355723 0.074003 0.318305 -0.326859 
2 0.987862 0.356105 0.074003 0.206074 -0.891025 
3 0.969875 0.365026 0.074003 0.133790 -0.930704 
4 0.944210 0.378836 0.074003 0.082276 -0.957979 
5 0.911495 0.394973 0.074003 0.043033 -0.978247 
6 0.872381 0.412926 0.074003 0.012034 -0.993965 
7 0.827561 0.432568 0.074003 -0.012724 -1.006342 
8 0.777784 0.453710 0.074003 -0.032414 -1.016078 
9 0.723850 0.476112 0.074003 -0.048010 -1.023724 
10 0.666609 0.500047 0.074003 -0.059905 -1.029517 
11 0.606945 0.525455 0.074003 -0.068405 -1.033637 
12 0.545769 0.552654 0.074003 -0.073563 -1.036129 


126 


: MISES 8.4% THICKNESS (COORDINATES SRE INPUT BY USER) 
: NLOWER = 25 


INCRESSae 


0.000000 


ley Of4emto Ome 0.074008 -0.075309 -1.036971 
Pe Oemeeeool)6= 6s 0.074008 -0.073531 -1.0356114 
15 0.362439 0.546480 0.074003 -0.067928 -1.033406 
16 0.304449 0.683297 0.074003 -0.057668 -1.028430 
mee Ovc2Geee> Ome > 0.074003 -0.041891 -1.029731 
Bey Ogbeseos Ot@eees 0.074003 -0.019114 -1.009512 
To) OF Mete2O O-etewe? 0.074803 0.015374 -0.993290 
20) Of WevevG O.e7eme2 0.074003 0.059669 -0.969707 
ile OND? Ofsoee? «8©60.074003 0.127774 -0.933931 
Gam Oeteeess Ltda > 8.074003 0.232984 -0.875795 
Jom Oedgezoo lols ©.074002 0.409236 -0.766612 
game OnOemece, ligegeem 0.074003 0.727179 -0.522323 
25 0.001884 1.653644 0.074003 0.945112 0.234282 
Gow O-O001ses2 0.545367 0.074003 -0.815326 1.347341 
27 0.009387 -0.275210 0.074003 -1.084184 1.443670 
Mo 0.024204 -0.497255 0.074003 -0.87Z601 1.368430 
fm © O2e255 -Ovoeter4 0.074003 -0.723499 1.312821 
Ee O20 525) -O.60e0e2 @.074003 -0.624167 1.274428 
om Ovimeo06 -0.62g609 0.074003 -0.552954 1.246176 
Be O-1ole4o -0.645575% 0.074003 -0.498371 1.224080 
eo ©.196e63 -0.649755 0.074003 -0.453794 1.205734 
Be ©.249455 -O.650079 0.074003 -0.415592 1.189787 
ps, 9.304240 -0.659575 ©.074003 -0.381557 1.175397 
Oa O.s024oq -OJee7ez. 0.0/4003 -0.34995/7 1.161877 
ey 0.422566 -0.64¢8020 0.074003 -0.319875 1.148858 
Soe O.4c7005 -O%eeeze. 0.074003 -0.290579 1.136037 
Bee Onoteyoo ~Oneeeee 2 © .0/74003 -0.267352 1.123099 
moe ©.605041 -0.642250 0.074003 -0.237604 1.109776 
pee O.cooo05 -0-641e95 0.@/4003 -0.200941 1.095875 
mee Ueigeenro “OnGe57G0 0.074003 -0.166763 1.081094 
pee O.//i¢o2 -Ovazooue 0.074003 -0.134640 1.065195 
Pee Ulso2 tol “Oey 0.074003 -0.0976/6 1.047701 
pe Onceezeo) -O7G2me) 0.074003 -0.056664 1.028039 
poe Ore o> -0.6547e0 0.074003 -0.010900 1.005436 
47 0.944210 -0.632196 0.074003 0.042413 0.978564 
48 0.969875 -0.629793 0.074003 0.107609 0.944665 
Poe, Usdereee -Oveeeeo0 0.074003 0.192062 0.898297 
ow 60..997429 -0.619634 0.074003 0.386207 0.826858 


oD =) 02000829 ep =" 0.203076 Gie= ~0.080325 


KEEKKKKKKKKKKRKKKKRKKRKRKRKRKRKRKRKAKRKRKRERERKKKRKRRK 


*xkX BEGIN UNSTEADY FLOW SOLUTION **%x 
KRAKKKKKAKARAKKAAKKARKRARARKRARARRRAKRAKA 


TIME STEP TK = 0.050000 TK - TKM1 = 0.050000 

RMPHA(T) = 2.516295 OMEGA(T) = -0.011248 

U(T) = 0.000000 ar) = -9.005624 

NITR VEW VYW WAKE THETA GAMMA 

0 0.999048 0.043619 0.050000 0.043633 0.740032E-01 
1 0.907832 0.005991 0.045393 0.006600 0.744799E-01 
2 0.904138 0.007297 0.045208 0.008070 0.744662E-01 
3 0.903985 0.007241 0.045201 0.008010 0.744652E-01 

CONVERGED SOLUTION OBTAINED AFTER NITR = 3 
J XG) Q(J) GAMMA crip v(J) 


1 0.997429 0.435837 0.074466 0.299470 -0.838202 
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2 0.987862 0.422852 0.074466 0.196120 -0.899569 
3 9.969875 0.421798 0.074466 0.132077 -0.936976 
4 9.344210 0.427256 0.074466 0.088439 -0.962416 
5 0.911495 0.435574 0.074466 0.055859 -0.981181 
3 9.872381 0.447181 0.074466 0.029822 -0.995676 
7 0.827561 0.460598 0.074466 0.008149 -1.007037 
3 0.777784 0.475913 0.074466 -0.010442 -1.015962 
9 0.723850 0.492830 0.074466 -0.026773 -1.022942 
10 9.866609 0.511571 0.074466 -0.041132 -1.028238 
i 9.606945 0.532044 0.074466 -0.053494 -1.032001 
-2 9.545769 0.554564 0.074466 -0.063531 -1.034262 
13 0.484008 0.579201 0.074466 -0.070990 -1.035018 
14 0.422591 0.606206 0.074466 -0.075185 -1.034204 
15 0.362439 0.635915 0.074466 -0.075487 -1.031672 
16 0.304449 0.669133 0.074466 -0.070722 -1.027008 
17 0.249485 0.706140 9.074466 -0.959745 -1.019773 
13 0.198363 0.748100 0.074466 -0.040826 -1.009171 
19 9.151840 0.796683 9.074466 -0.011154 -0.993762 
20 0.110606 9.853825 0.074466 0.033396 -0.971231 
21 0.075269 0.924099 0.074466 0.100715 -0.936853 
22 0.046353 1.014820 0.074466 0.205876 -0.880676 
23 0.024285 1.143358 0.074466 0.382654 -0.776542 
24 0.009388 1.351849 0.974466 0.703606 -0.535871 
25 0.001884 1.634712 0.074466 0.952740 0.209596 
26 0.001884 0.564272 0.074466 -0.746972 1.322641 
27 0.909387 -0.246433 0.074466 -1.037466 1.430099 
28 0.024284 -0.4678138 0.074466 -0.838061 1.360476 
28 0.046353 -0.551795 0.074466 -0.693563 1.307915 
30 0.075269 -0.590860 0.074466 -0.596473 1.271481 
31 0.110606 -0.611392 0.074466 -0.527117 1.244626 
32 9.151840 -0.622549 0.074466 -0.474818 1.223580 
33 0.198363 -0.629333 0.074466 -0.433321 1.206047 
34 0.249485 -0.633515 0.074466 -0.399064 1.190716 
33 0.304448 -0.636433 0.074466 -0.369826 1.176793 
36 0.362436 -0.539059 0.074466 -0.343641 1.163586 
37 0.422588 -0.641343 0.074466 -0.319346 1.150745 
38 0.484005 -0.643768 0.074466 -0.295850 1.137963 
39 0.545766 -0.646483 0.074466 -0.272088 1.124936 
40 0.606941 -0.649578 0.074466 -0.247068 1.111386 
41 0.666606 -0.652918 0.074466 -0.220048 1.097122 
42 0.723848 -0.656699 0.074466 -0.190134 1.081849 
43 0.777782 -0.660707 0.074466 -0.156552 1.065285 
44 0.827561 -0.665236 0.074466 -0.118313 1.046980 
45 0.872381 -0.670135 0.074466 -0.074279 1.026307 
46 0.911495 -0.675261 0.074466 -0.023233 1.002476 
47 0.944210 -0.680614 0.074466 0.036802 0.974108 
48 0.969875 -0.686543 0.074466 0.109850 0.938386 
49 0.987862 -0.692719 0.074466 0.203356 0.889792 
50 0.997429 -0.699982 0.074466 0.333160 0.815602 
CD = 0.001539 CL = 0.302054 CM = -0.088450 
TRAILING VORTICES DATA 
M = -K(M) ¥(M) CIRC 
1 1.022599 9.000131 -9.000933 
TIME STEP TK = 0.749999 TK - TKM1 = 0.050000 
ALPHA(T) = 4.999996 OMEGA(T) = -0.087266 
U(T) = 0.000000 V(T) = -0.043633 
NITR VW VYW WAKE THETA GAMMA 


O 0.905684 =0.000916 0.045284 -0,0010RZ “O-2Ge2eSEr00 
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1 0.905735 -0.000649 0.045287 -0.000717 
CONVERGED SOLUTION OBTAINED AFTER NITR= 1 

J ea) Q(3) GAMMA CP(a) 

1 0.997429 1.115997 0.106565 0.311649 

2 0.987862 1.060333 0.106565 0.221106 

3 0.969875 1.031653 0.106565 0.170306 

4 0.944210 1.013085 0.106565 0.141163 

5 0.911495 0.998141 0.106565 0.123949 

6 0.872381 0.985265 0.106565 0.114073 

7 0.827561 0.974241 0.106565 0.109016 

8 0.777784 0.964993 0.106565 0.107136 

9 0.723850 0.957451 0.106565 0.107130 
10 0.666609 0.952120 0.106565 0.108458 
11 0.606945 0.949235 0.106565 0.110657 
12 0.545769 0.949393 0.106565 0.113665 
13 0.484008 0.952960 0.106565 0.117540 
14 0.422591 0.960451 0.106565 0.122480 
15 0.362439 0.972455 0.106565 0.128967 
16 0.304449 0.990059 0.106565 0.138078 
17 0.249485 1.013807 0.106565 0.150962 
18 0.198363 1.045081 0.106565 0.169616 
19 0.151840 1.085785 0.106565 0.197364 
20 0.110606 1.133024 0.106565 0.239070 
21 0.075269 1.206453 0.106565 0.303822 
22 0.046353 1.298127 0.106565 0.407818 
23 0.024285 1.428849 0.106565 0.583190 
24 0.009388 1.631804 0.106565 0.872165 
25 0.001884 1.819807 0.106565 0.765241 
26 0.001884 0.373179 0.106565 -1.559307 
27 0.009387 -0.529374 0.106565 -1.564167 
28 0.024284 -0.755145 0.106565 -1.202155 
29 0.046353 -0.836350 0.106565 -0.991312 
30 0.075269 -0.874103 0.106565 -0.864650 
31 0.110606 -0.896239 0.106565 -0.781035 
32 0.151840 -0.912096 0.106565 -0.721021 
33 0.198363 -0.926607 0.106565 -0.674005 
34 0.249485 -0.941343 0.106565 -0.634257 
35 0.304448 -0.957406 0.106565 -0.598321 
36 0.362436 -0.975544 0.106565 -0.563539 
37 0.422588 -0.995441 0.106565 -0.528592 
38 0.484005 -1.017302 0.106565 -0.492391 
39° 0.545766 -1.041010 0.106565 -0.454043 
40 0.606941 -1.066387 0.106565 -0.412924 
41 0.666606 -1.093023 0.106565 -0.368740 
42 0.723848 -1.120818 0.106565 -0.321026 
43 0.777782 -1.149241 0.106565 -0.269645 
44 0.827561 -1.178310 0.106565 -0.213999 
45 0.872381 -1.207644 0.106565 -0.153563 
46 0.911495 -1.236895 0.106565 -0.087684 
47 0.944210 -1.265976 0.106565 -0.014749 
48 0.969875 -1.296105 0.106565 0.069103 
49 0.987862 -1.330154 0.106565 0.171240 
50 0.997429 -1.380203 0.106565 0.308161 

Chw= Ge0e9'687 CL = 0.645338 CM = 

TRAILING VORTICES DATA 

M X(M) Y(M) CIRC 

1 1.700760 0.077052 -0.000933 

2 1.651289 0.072156 -0.001625 

3 1.602002 0.066331 -0.002229 

4 1.552845 0.060031 -0.002784 

5 1.503779 0.053471 -0.003304 

6 1.454797 0.046788 -0.003797 

7 1.405895 0.040107 -0.004256 


IZ 


DODD OKRPBBBBHBHBHBPBHBHBHHBPHPBPHPHHHO 


.106563E+00 


oe 


Re OrsyWans, 
~Jomes3 
sJOS soo 
soo /6 
.008098 
eOls416 
-016142 
OOS 
.016594 
015062 
,OLZoe2 
JOSS 
.004971 
sao oa0D 
~992654 
383355 
~972486 
.957451 
936848 
e207 oS 
-863894 
./ 25034 
isis 2 
- 369636 
-485826 
<a lO26 
2908s! 
.468068 
wCooee 
~ 338450 
30214 
~274529 
-251843 
(Zo, 
-214148 
. 35805 
eed 
~162518 
» 144540 
Zo 2 
Ooo 
083543 
OSes 
~034027 
wo Zo9 
soc l00 
Jc 
Joong 
Osage! 
- 746084 


8 1.357057 0.033554 -0.004687 
9 1.308305 0.027216 -0.005100 
10 1.259669 0.021159 -0.005481 
11 1.211178 0.015541 -0.005801 
12 1.162916 0.010481 -0.006100 
13. 1.115048 0.006079 -0.006349 
14 1.067918 0.002406 -0.006559 
15 1.022643 -0.000016 -0.006720 
TIME STEP TK = 1.449992 TK - TKM1 = 0.050000 
ALPHA(T) = 7.483698 OMEGA(T) = -0.011249 
uU(T) = 0.000000 V(T) = -0.005625 
NITR VXW VYW WAKE THETA GAMMA 
0 0.901699 0.009872 0.045088 0.010948 0.145377E+00 
1 0.901200 0.011028 0.045063 0.012236 0.146997E+00 
CONVERGED SOLUTION OBTAINED AFTER NITR = 1 
J K(J) Q(J) GAMMA CP(J) Vien 
1 0.997429 1.101898 0.146996 0.332408 -0.864078 
2 0.987862 1.099609 0.146996 0.228004 -0.921138 
3 0.969875 1.116333 0.146996 0.160643 -0.955077 
4 0.944210 1.140985 0.146996 0.114868 -0.976579 
5 0.911495 1.168619 0.146996 0.082627 -0.990894 
6 0.872381 1.197904 0.146996 0.060528 -1.000257 
7 0.827561 1.228671 0.146996: 0.046710 -1.005878 
8 0.777784 1.260826 0.146996 0.039936 -1.008490 
9 0.723850 1.294203 0.146996 0.039112 -1.008660 
10 0.666609 1.329222 0.146996 0.043766 -1.006560 
11 0.606945 1.366027 0.146996 0.053332 -1.002352 
12 0.545769 1.405159 0.146996 0.067653 -0.995947 
13 0.484008 1.446882 0.146996 0.086578 -0.987214 
14 0.422591 1.491619 0.146996 0.110103 -0.975912 
15 0.362439 1.539856 0.146996 0.138557 -0.961606 
16 0.304449 1.592619 0.146996 0.172964 -0.943442 
17 0.249485 1.650390 0.146996 0.214399 -0.920465 
18 0.198363 1.714437 0.146996 0.265094 -0.890942 
19 0.151840 1.786607 0.146996 0.328746 -0.851951 
20 0.110606 1.868882 0.146996 0.410545 -0.798852 
21 0.075269 1.965495 0.146996 0.519597 -0.722328 
22 0.046353 2.082446 0.146996 0.668455 -0.603478 
23 0.024285 2.230950 0.146996 0.866926 -0.394782 
24 0.009388 2.419835 0.146996 1.009466 0.050859 
25 0.001884 2.339782 0.146996 -0.471799 1.214887 
26 0.001884 -0.156346 0.146996 -4.389847 2.320095 
27 0.009387 -1.322127 0.146996 -3.033496 2.003043 
28 0.024284 -1.560176 0.146996 -2.015200 1.727204 
29. 0.046353 -1.622657 0.146996 -1.505832 1.569752 
30 0.075269 -1.634560 0.146996 -1.212790 1.470549 
271 0.110606 -1.528102 0.146996 -1.021804 1.401554 
32 9.151340 -1.513625 0.146996 -0.3856i6 1.550007 
33. 0.198363 -1.596423 0.146996 -0.730655 i.305004 
34 0.249485 -1.578180 0.146996 -0.695048 1.274895 
35 0.304448 -1.560042 0.146996 -0.621833 1.245407 
36 0.362436 -1.542861 0.146996 -0.556431 1.218916 
37 0.422588 -1.526391 0.146996 -0.496616 1.194569 
38 0.484005 -1.510876 0.146996 -0.440592 1.171595 
39 0.545766 -1.496317 0.146996 -0.387199 1.149437 
40 0.606941 -1.482638 0.146996 -0.335611 1.127617 
41 0.666606 -1.469497 0.146996 -0.285518 1.105863 
42 0.723848 -1.456874 0.146996 -0.236273 1.083745 
43 0.777782 -1.444335 0.146996 -0.187542 1.060991 


44 0.827561 -1.431955 0.146996 -0.138458 1.037087 
45 0.372381 -1.419472 0.146996 -0.088061 1.011468 
46 0.911495 -1.406547 0.146996 -0.035055 0.9834i1 
47 0.944210 -1.393024 0.146996 0.023027 0.951546 
48 0.969875 -1.379868 0.146996 0.091342 0.912885 
49 0.987862 -1.368515 0.146996 0.178655 0.861777 
50 0.997429 -1.365097 0.146996 0.303827 0.784388 
CD = 0.030956 ci 0.713821 CM = -0.190685 
TRAILING VORTICES DATA 
M = X(M)- ¥(M) CIRC 
1 2.383780 0.232008 -0.000933 
2 2.333846 0.225080 -0.001625 
3 2.284404 0.216393 -0.002229 
4 2.235273 0.206695 -0.002784 
5 2.186347 0.196301 -0.003304 
6 2.137600 0.185376 -0.003797 
7 2.089004 0.174067 -0.004256 
8 2.040442 0.162550 -0.004687 
9 1.991957 0.150853 -0.005100 
10 1.943572 0.138989 -0.005481 
11 1.895155 0.127193 -0.605801 
12 1.846638 0.115584 -0.006100 
13. 1.798126 0.104170 -0.006349 
14 1.749496 0.093078 -0.006559 
15 1.700793 0.082359 -0.006720 
16 1.651988 0.072087 -0.006839 
17 1.603077 0.062348 -0.006896 
18 1.554079 0.053196 -0.006916 
19 1.505019 0.044648 -0.006885 
20 1.455917 0.036758 -0.006795 
21 1.406791 0.029591 -0.006642 
22 1.357685 0.023172 -0.006443 
23 1.308651 0.017529 -0.006177 
24 1.259734 0.012684 -0.005852 
25 1.211021 0.008642 -0.005465 
26 1.162620 0.005401 -0.005014 
27 1.114706 0.002948 -0.004498 
28 1.067624 0.001210 -0.003917 
29° 1.022530 0.000276 -0.003269 
TIME STEP TK = 1.999990 TK - TKM1 = 0.200000 
ALPHA(T) = 7.500000 OMEGA(T) = 0.000000 
U(T) = 0.000000 vV(T) = 0.000000 
NITR VXW VYW WAKE THETA GAMMA 
O 0.939783 0.026143 0.188029 0.027811 0.156187E+00 
1 0.948367 0.030228 0.139770 0.031863 0.160749E+00 
2 0.948647 0.030081 0.189825 0.031699 0.160765£+00 
CONVERGED SOLUTION OBTAINED AFTER NITR = 2 
J ap Q(J) GAMMA CP(J) ese 
1 0.997429 1.116135 0.160766 0.320772 -0.851401 
2 0.987862 1.110748 0.160766 0.221351 -0.907726 
3 0.969875 1.126034 0.160766 0.158555 -0.941336 
4 0.944210 1.150497 0.160766 0.116479 -0.962935 
5 0.911495 1.179177 0.160766 0.086724 -0.977640 
6 0.872381 1.210700 0.160766 0.065709 -0.987588 
7 0.827561 1.244766 0.160766 0.051593 -0.993866 
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